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ABSTRACT 

Malyshkina, Nataliya V. Ph.D., Purdue University, December 2008. Markov Switch- 
ing Models: an Apphcation to Roadway Safety. Major Professors: Fred L. Mannering 
and Andrew P. Tarko. 

In this research, two-state Markov switching models are proposed to study accident 
frequencies and severities. These models assume that there are two unobserved states 
of roadway safety, and that roadway entities (e.g., roadway segments) can switch 
between these states over time. The states are distinct, in the sense that in the 
different states accident frequencies or severities are generated by separate processes 
(e.g., Poisson, negative binomial, multinomial logit). Bayesian inference methods and 
Markov Chain Monte Carlo (MCMC) simulations are used for estimation of Markov 
switching models. To demonstrate the applicability of the approach, we conduct the 
following three studies. 

In the first study, two-state Markov switching count data models are considered as 
an alternative to zero-inflated models, in order to account for preponderance of zeros 
typically observed in accident frequency data. In this study, one of the states of road- 
way safety is a zero-accident state, which is perfectly safe. The other state is an un- 
safe state, in which accident frequencies can be positive and are generated by a given 
counting process - a Poisson or a negative binomial. Two-state Markov switching 
Poisson model, two-state Markov switching negative binomial model, and standard 
zero-inflated models are estimated for annual accident frequencies on selected Indiana 
interstate highway segments over a five-year time period. An important advantage of 
Markov switching models over zero-inflated models is that the former allow a direct 
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statistical estimation of what states specific roadway segments are in, wfiile the later 
do not. 

In the second study, two-state Markov switching Poisson model and two-state 
Markov switching negative binomial model are estimated using weekly accident fre- 
quencies on selected Indiana interstate highway segments over a five-year time period. 
In this study, both states of roadway safety are unsafe. Thus, accident frequencies 
can be positive and are generated by either Poisson or negative binomial processes in 
both states. It is found that the more frequent state is safer and it is correlated with 
better weather conditions. The less frequent state is found to be less safe and to be 
correlated with adverse weather conditions. 

In the third study, two-state Markov switching multinomial logit models are esti- 
mated for severity outcomes of accidents occurring on Indiana roads over a four-year 
time period. It is again found that the more frequent state of roadway safety is corre- 
lated with better weather conditions. The less frequent state is found to be correlated 
with adverse weather conditions. 

One of the most important results found in each of the three studies, is that in 
each case the estimated Markov switching models are strongly favored by accident 
frequency and severity data and result in a superior statistical fit, as compared to the 
corresponding standard (single-state) models. 
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CHAPTER 1. INTRODUCTION 



This chapter explains the motivation and objectives of the present research, and the 
organization of this dissertation. 



1.1 Motivation and research objectives 



According to Bureau of Transportation Statistics [BTS|, |2008| . in 2006, 99.55% 
of all transportation related accidents (including air, railroad, transit, waterborne 
and pipeline accidents) were motor vehicle accidents on roadways. Motor vehicle 
accidents result in fatalities, injuries and property damage, and represent high cost 
not only for involved individuals but also for our society as a whole. In particular, 
on average, about one-quarter of the costs of crashes is paid directly by the party 
involved, while the society pays the rest. As an example of the economic burden 
related to motor vehicle crashes, in the year 2000 the estimated cost of accidents 
occurred in the United States was 231 billion doll ars, which is about 820 dollars per 



BTS 



20081]. These numbers show 



person or 2 percent of the gross domestic product 
that roadway vehicle travel safety has an enormous importance for our society and for 
the national economy. As a result, extensive research on roadway safety is ongoing, 
in order to better understand the most important factors that contribute to vehicle 
accidents. 

In general, there are two measures of roadway safety that are commonly consid- 
ered: 



1. The first measure evaluates accident frequencies on roadway segments. Accident 
frequency on a roadway segment is obtained by counting the number of acci- 
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dents occurring on this segment during a specified period of time. Then count 
data statistical models (e.g. Poisson, negative binomial models and their zero- 
inflated counterparts) arc estimated for accident frequencies on different road- 
way segments. The explanatory variables used in these models are the roadway 
segment characteristics (e.g. roadway segment length, curvature, slope, type, 
pavement quality, etc). 

2. The second measure evaluates accident severity outcomes as determined by the 
injury level sustained by the most severely injured individual (if any) involved 
into the accident. This evaluation is done by using data on individual accidents 
and estimating discrete outcome statistical models (e.g. ordered probit and 
multinomial logit models) for the accident severity outcomes. The explanatory 
variables used in these models are the individual accident characteristics (e.g. 
time and location of an accident, weather conditions and roadway characteristics 
at the accident location, characteristics of the vehicles and drivers involved, etc). 

These two measures of roadway safety are complementary. On one hand, an 
accident frequency study provides a statistical model of the probability of an accident 
occurring on a roadway segment. On the other hand, an accident severity study 
provides a statistical model of the conditional probability of a severity outcome of an 
accident, given the accident occurred. The unconditional probability of the accident 
severity outcome is the product of its conditional probability and the probability of 
the accident. 

The main objective of this research study is to propose a new statistical approach 
to modeling accident frequencies and severities, which may provide new guidance to 
theorists and practitioners in the area of roadway safety. Our approach is based on 
application of two-state Markov switching models of accident frequencies and severi- 
ties. These models assume an existence of two unobserved states of roadway safety. 
The roadway entities (e.g., roadway segments) arc assumed to be able to switch be- 
tween these states over time, and the switching process is assumed to be Markovian. 
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Accident frequencies and severity outcomes are assumed to be generated by two dis- 
tinct data-generating processes in the two states. Two-state Markov switching models 
avoid several drawbacks of the popular conventional models of accident frequencies 
and severities. We estimate Markov switching models and compare them to the con- 
ventional models. We find that the former are strongly favored by accident frequency 
and severity data and provide a superior statistical fit as compared to the later. Be- 
cause of the complexity of Markov switching models, this research employs Bayesian 
inference and Markov Chain Monte Carlo (MCMC) simulations for their statistical 
estimation. 



1.2 Organization 



An overview of the previous research on accident frequency and severity is pre- 
sented in Chapter [21 Chapter [3] gives specification of the two-state Markov switching 
and conventional models that are proposed, considered and estimated in this study. 
Bayesian inference methods are given in Chapter HI Chapter [5] presents Markov Chain 
Monte Carlo (MCMC) simulation techniques used for Bayesian inference and model 
estimation in this study. The model estimation results for accident frequencies are 
presented in Chapter [61 The model estimation results for accident severities are given 
in Chapter [71 Finally, we discuss our results and give conclusions in Chapter [HI Some 
of the results are given in the Appendix at the end. 
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CHAPTER 2. LITERATURE REVIEW 



This chapter includes a brief overview of the previous roadway safety studies of ac- 
cident frequencies and severities. First, we give an overview of accident frequency 
studies and standard statistical models used for accident frequencies. Then we re- 
view previous work on severities of accidents. Finally, we discuss studies that consider 
both accident frequencies and accident severities. The literature review of this chap- 
ter does not claim to be full or exhaustive. A more detailed literature review, as 
well as a comprehensive description of c onventional methodologie s commonly used in 
roadway safety studies, can be found in [Washington et al.l [2003|. 



2.1 Accident frequency studies 

Considerable research has been conducted on understanding and predicting ac- 
cident frequencies (the number of accidents occurring on roadway segments over a 
given time period). Because accident frequencies are non-negative integers, count 
data models are a reasonable statistical modeling approach. Simple modeling ap- 
proaches include Poisson models and negative binomial (NB) models. These models 
assume a single process for accident data generation (a Poisson process or a negative 
binomial process) and involve a nonlinear regression of the observed accident fre- 
quencies on various roadway-segment characteristics (such as roadway geometric and 
environmental factors). Selected previous research on accident frequencies, conducted 
by application of count data models, is as follows: 



Hadi et al. 



1995l | used negative binomial models to estimate the effect of cross 



section roadway design elements (e.g. presence of curb, lane width) and traffic 
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volume on accident frequencies for different types of highways. The authors 
found that some cross section design elements can influence accident rates (e.g. 
lane width, interchange presence, speed limit) and that some other do not have 
any effect on number of accidents (e.g. type of friction course material). 



Shankar et al 



1995l | applied a negative binomial model to an accident data col- 
lected in Washington State. Roadway geometries of fixed-equal length roadway 
segments (e.g. horizontal and vertical alignments), weather, and other seasonal 
effects were analyzed along with overall accident frequencies of specific accident 
types (e.g., rear-end and same direction accidents). This research concluded 
that highway segments with challenging geometries as well that fre- 

quently experience adverse weather conditions are important determinants of 
accident frequency. 



Poch and ManneringI 1996l | estimated a negative binomial regression of the fre- 
quencies of accidents at intersection approaches in Seattle suburban areas. The 
authors of this paper considered traffic volume, geometric characteristics of in- 
tersection approaches (e.g. approach sight-distance, speed limit) and approach 
signalization characteristics (e.g. eight-phase signal) as the model explanatory 
variables. Authors found a significant influence of some of these variables on 
accident frequencies at intersection approaches. In particular, they found that 
high left-turn and opposite traffic volumes considerably increase numbers of 
accidents at intersection approaches. 



Miaou and Lord 



20031], based on accident data collected in Toronto, examined 
generally accepted statistical models (Poisson and NB) applied to accident fre- 
quencies intersections. By using the empirical Bayes method, mathematical 
properties and performance of different popular model functional forms were 
considered. The authors questioned invariability of the dispersion parameter, 
given the complexity of the traffic interaction in an intersection area. In addi- 
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tion, the full Bayes statistical approach was used for model specification and 
estimation. 



Park and LordI 2008[ | recently considered finite mixture Poisson and negative 
binomial models of accident frequencies, in order to account for heterogenous 
populations of accident data. Accident data heterogeneity can result from 
data generation by distinct (Poisson or N B) processes that oper ate in differ- 



Park and LordI |2008| suggested a 



ent unobserved states of roadway safety, 
two-component finite mixture negative binomial model as the best model to 
account for accident data heterogeneity in their data sample. 



Recently, lAnastasopoulos and ManneringI |2008l | applied random parameters 
count models to the analysis of accident frequencies. The authors found these 
models to be beneficial for accident frequency prediction. Random parame- 
ter models can potentially define unique parameters for each roadway segment, 
but these models still assume a single state for each segment. This single- 
state assumption wo uld also be true for count models with random effects [see 



Shankar et al. 



1998 |. 



Anastasopoulos et al. 



2008| were the first to use tobit regression models for 
prediction of accident rates (accident rates are number of accidents happened 
per unit roadway segment length and per unit averaged annual daily traffic 
volume). They considered five-year accident data and found that international 
roughness index (of the pavement), pavement rutting, the pavement's condition 
rating, median types and width, shoulder widths, number of ramps and bridges, 
horizontal and vertical curves, rumble strips, annual average daily travel and 
the percent of combination truck in the traffic stream have a significant impact 
on accident rates. 

Because a preponderance of zero-accident observations is often observed in empir- 
ical data, some researchers have applied zero-inflated Poisson (ZIP) and zero-inflated 
negative binomial (ZINB) models for predicting accident frequencies. Zero-inflated 
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models assume a two-state process for accident data generation. One state is assumed 
to be perfectly safe with zero accidents (over the duration of time being considered). 
The other state is assumed to be unsafe with a possibility of nonzero accident frequen- 
cies in which accidents can happen and accident frequencies are generated by some 
given counting process (Poisson or negative binomial). Below are selected studies 
that are based on an application of zero-inflated count data models: 



Miaoul |1994| applied Poisson regression, zero-inflated Poisson (ZIP) regression, 
and NB regression to determine a relationship between geometric design char- 
acteristics of roadway segments and the number of truck accidents. Results 
suggest that under the maximum likelihood estimation (MLE) method, all 
three models perform similarly in terms of estimated truck-involved accident 
frequencies across roadway segments. To model the relationship, the author 
recommended the use of a Poisson regression as an initial model, then the use 
of a negative binomial model if the accident frequency data is overdispersed, 
and the use of a zero-inflated Poisson model if the data contains an excess of 
zero observations. 



Shankar et al. 



1997( 1 studied the distinction between safe and unsafe roadway 
segments by estimating zero-inflated Poisson and zero-inflated negative binomial 
models for accident frequencies in Washington State. The authors established 
the underlying principles of zero-inflated models, based on a two-state data- 
generating process for accident frequencies. The two states are a safe state that 
corresponds to the zero accident likelihood on a roadway segment, and an unsafe 
state. The results show that two-state zero-inflated structure models provide 
a superior statistical fit to accident frequency data as compared to the con- 
ventional single-state models (without zero- inflation). Thus, the authors found 
that zero-inflated models are helpful in revealing and understanding important 
factors that affect accident frequencies with preponderance of zeros. 
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Lord et al.l |2005l . l2007l | addressed the question of choosing the best approach 
to the modehng of roadway accident data by using count data models (e.g. 
whether to use standard single-state or zero-inflated models). Authors argued 
that an application of zero-inflated models to the analysis of accident data with 
a preponderance of zeros is not a defensible modeling approach. They argued 
that an excess of zeros can be caused by an inappropriate data collection and 
by many other factors, instead of due to a two-states process. In addition, 
they claimed that it is unreasonable to expect some roadway segments to be 
always perfectly safe and questioned "safe" and "unsafe" state definitions. The 
authors also argued that zero-inflated models do not explicitly account for a 
likely po ssibility for roadway se gments to change in time from one state to 



another. 



Lord et al. 



2005 



20071 ] concluded that, while an application of zero- 



inflated models often provides a better statistical fit to an observed accident 
frequency data, the applicability of these models can be questioned. 



2.2 Accident severity studies 



Research efforts in predicting accident severity, such as property damage, injury 
and fatality, are clearly very important. In the past there has been a large number 
of studies that focused on modeling accident severity outcomes. The probabilities of 
severity outcomes of an accident are conditioned on the occurrence of the accident. 
Common modeling approaches of accident severity include multinomial logit models, 
nested logit models, mixed logit models and ordered probit models. All accident 
severity models involve nonlinear regression of the observed accident severity out- 
comes on various accident characteristics and related factors (such as roadway and 
driver characteristics, environmental factors, etc). Some of the past accident severity 
studies follows: 
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O'Donnell and Connoii 1996| explored severity of motor vehicle accidents in 
Australia by estimating the parameters of ordered multiple choice models: or- 
dered logit and probit models. By studying driver, passengers and vehicle char- 
acteristics (e.g. vehicle type, seating position of vehicle occupants, blood alcohol 
level of a driver), the authors found the effects of these characteristics on the 
probabilities of different types of severity outcomes. For example, they found 
that the older the victims are and the higher the vehicle speeds are, the higher 
the probabilities of serious injuries and deaths are. 



Shankar and Mannering] 1996| estimated the likelihoods of motorcycle rider ac- 
cident severity outcomes. In their research work, a multinomial logit model was 
applied to a 5-year Washington state data for single-vehicle motorcycle colli- 
sions. It was found that a helmeted-riding is an effective means of reducing 
injury severity in any types of collisions, except in fixed-object collisions. At 
the same time, alcohol-impaired riding, high age of a motorcycle rider, ejection 
of a rider, wet pavement, interstate as a roadway type, speeding and rider inat- 
tention were found to be the factors that increase roadway motorcycle accident 
severity. 



Shankar et al.l 1996| used a nested logit model for statistical analysis of acci- 
dent severity outcomes on rural highways in Washington State. They found that 
environment conditions, highway design, accident type, driver and vehicle char- 
acteristics significantly influence accident severity. They found that overturn 
accidents, rear-end accidents on wet pavement, fixed-object accidents, and fail- 
ure to use the restraint belt system lead to higher probabilities of injury or /and 
fatality accident outcomes, while icy pavement and single-vehicle collisions lead 
to higher probability of property damage only outcomes. 



Duncan et al. 



1998l | applied an ordered probit model to injury severity out- 



comes in truck-passenger car rear-end collisions in North Carolina. They found 
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that injury severity is increased by darkness, high speed differentials, high speed 
hmits, wet grades, drunk driving, and being female. 



Chang and Mannering] 1999l | focused on the effects of trucks and vehicle occu- 
pancies on accident severities. They estimated nested logit models for sever- 
ity outcomes of truck-involved and non-truck-involved accidents in Washington 
State and found that accident injury severity is noticeably worsened if the ac- 
cident has a truck involved, and that the effects of trucks are more significant 
for multi-occupant vehicles than for single-occupant vehicles. 



Khattak 



2001| estimated ordered probit models for severity outcomes of multi- 



vehicle rear-end accidents in North Carolina. In particular, the results of his 
research indicate that in two- vehicle collisions the leading driver is more likely 
to be severely injured, in three- vehicle collisions the driver in the middle is more 
likely to be severely injured, and being in a newer vehicle protects the driver in 
rear-end collisions. 



UlfarssonI 200 1| . lUlfarsson and Mannering! 2004l | focused on male and female 
differences for accident severity outcomes. They used multinomial logit models 
and accident data from Washington State. They found significant behavioral 
and physiological differences between genders, and also found that probability 
of fatal and disabling injuries is higher for females as compared to males. 



Kockelman and KweonI 2002| applied ordered probit models to modeling of 
driver injury severity outcomes. They used a nationwide accident data sample 
and found that pickups and sport utility vehicles are less (more) safe than 
passenger cars in single-vehicle (two-vehicle) collisions. 



Khattak et al. 



2002l | focused on the safety of aged drivers in the United States. 
Nine-year Iowa-statewide accident data was considered and the ordered probit 
modeling technique was implemented for accident severity modeling. Authors 
inspected vehicle, roadway, driver, collision, and environmental characteristics 
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as factors that may potentially effect accident severity of aged drivers. The 
modeling results were consistent with a common sense, for example, an animal- 
related accident tends to have severe consequences for elderly drivers. Also, 
it was found that accidents with farm vehicles involved are highly severe for 
elderly drivers in Iowa. 



Abdel-Atyl 2003| used ordered probit models for analysis of driver injury sever- 
ity outcomes at different road locations (roadway segments, signalized intersec- 
tions, toll plazas) in Central Florida. He found higher probabilities of severe 
accident outcomes for older drivers, male drivers, those not wearing seat belt, 
drivers who speed, those who drove vehicles struck at the driver's side, those 
who drive in rural areas, and drivers using electronic toll collection device (E- 
Pass) at toll plazas. 



Yamamoto and Shankarl 2004l | applied bivariate ordered probit models to an 
analysis of driver's and passenger's injury severities in collisions with fixed ob- 
jects. They considered a 4-year accident data sample from Washington State 
and found that collisions with leading ends of guardrail and trees tend to cause 
more severe injuries, while collisions with sign posts, faces of guardrail, concrete 
barrier or bridge and fences tend to cause less severe injuries. They also found 
that proper use of vehicle restraint system strongly decreases the probability of 
severe injuries and fatalities. 



Khorashadi et al. 



2005| explored the differences of driver injury severities in 
rural and urban accidents involving large trucks. Using four years of California 
accident data and multinomial logit model approach, they found considerable 
differences between rural and urban accident injury severities. In particular, 
they found that the probability of severe/fatal injury increases by 26% in rural 
areas and by 700% in urban areas when a tractor-trailer combination is involved, 
as opposed to a single-unit truck being involved. They also found that in ac- 
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ci dents where alcohol or drug use is identified, the probability of severe/fatal 
injury is increased by 250% and 800% in rural and urban areas respectively. 



Islam and Mannering 2006| studied driver aging and its effect on male and 
female single- vehicle accident injuries in Indiana. They employed multinomial 
logit models and found significant differences between different genders and 
age groups. Specifically, they found an increase in probabilities of fatality for 
young and middle-aged male drivers when they have passengers, an increase in 
probabilities of injury for middle-aged female drivers in vehicles 6 years old or 
older, and an increase in fatality probabilities for males older than 65 years old. 



Malyshkinal 2006| , iMalyshkina and Mannering 2006| focused on the relation- 
ship between speed limits and roadway safety. Their research explored the 
influence of the posted speed limit on the causation and severity of accidents. 
Multinomial logit statistical models were estimated for causation and severity 
outcomes of different types of accidents on different road classes. The results 
showed that speed limits do not have a statistically significant adverse effect 
on unsafe-speed-related causation of accidents on all roads. At the same time 
higher speed limits generally increase the severity of accidents on the majority 
of roads other than interstate highways (on interstates speed limits were found 
to have statistically insignificant effect on accident severity). 



Savolainen 



20061 ]. ISavolainen and Manneringl 2007l | focused on the important 



topic of motorcycle safety on Indiana roads. They used multinomial and nested 
logit models and found that poor visibility, unsafe speed, alcohol use, not wear- 
ing a helmet, right-angle and head-on collisions, and collisions with fixed objects 
increase severity of motorcycle-involved accidents. 



Milton et al. 



20081], by using accident severity data from Washington State, 
estimated a mixed logit model with random parameters. This approach allows 
estimated model parameters to vary randomly across roadway segments to ac- 
count for unobserved effects that can be related to other factors influencing 
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roadway safety. Authors found that, on one hand, some roadway characteristic 
parameters (e.g. pavement friction, number of horizontal curves) can be taken 
as fixed. On the other hand, other model parameters, such as weather effects 
and volume-related model parameters (e.g. truck percentage, average annual 
snowfall), are random and normally- distributed. 



Eluru and Bhat 



2007l | modeled a seat belt use endogeneity to accident severity 
due to unsafe driving habits of drivers not using seat belts. For severity out- 
comes, the authors considered a system of two mixed probit models with random 
coefficients estimated jointly for seat belt use dummy and severity outcomes. 
The probit models included random variables that moderate the influence of the 
primary explanatory attributes associated with drivers. The estimation results 
highlight the importance of moderation effects, seat belt use endogeneity and 
the relation of between failure to use seat belt and unsafe driving habits. 



2.3 Mixed studies 

Several previous research studies considered modeling of both accident frequen- 
cies and accident severity outcomes. It is beneficial to look at both frequencies and 
severities simultaneously because, as mentioned above, an unconditional probability 
of the accident severity outcome is the product of its conditional probability and the 
accident probability. Several mixed studies, which consider both accident frequency 
and severity, are as follows. 



Carson and Manneringj |2001| studied the effect of ice warning signs on ice- 
accident frequencies and severities in Washington State. They modeled accident 
frequencies and severities by using zero-inflated negative binomial and logit 
models respectively. They found that the presence of ice warning signs was not 
a signiflcant factor in reducing ice-accident frequencies and severities. 
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Lee and Mannering] [2002l | estimated zero-inflated count-data models and nested 
logit models for frequencies and severities of run-off-roadway accidents in Wash- 
ington State. They found that run-off-roadway accident frequencies can be re- 
duced by avoiding cut side slopes, decreasing (increasing) the distance from 
outside shoulder edge to guardrail (light poles), and decreasing the number 
of isolated trees along roadway. The results of their research also show that 
run-off-roadway accident severity is increased by alcohol impaired driving, high 
speeds, and the presence of a guardrail. 



Kweon and Kockelman 



2003l | studied probabilities of accidents and accident 
severity outcomes for a given fixed driver exposure (defined as the total miles 
driven). They used Poisson and ordered probit models, and considered a na- 
tionwide accident data sample. After normalization of accident rates by driver 
exposure, the results of their study indicated that young drivers are far more 
crash prone than other drivers, and that sport utility vehicles and pickups are 
more likely to be involved into rollover accidents. 
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CHAPTER 3. MODEL SPECIFICATION 



In this chapter we specify the statistical models that are used and estimated in the 
present study. First, we consider standard (conventional) models commonly used 
in accident studies. These are count data models for accident frequencies (Poisson, 
negative binomial models and their zero-inflated counterparts) and discrete outcome 
models for accident severity outcomes (multinomial logit models). Then we explain 
Markov process for the state of roadway safety. Finally, we present two-state Markov 
switching models for accident frequencies and severities. In each of the two states the 
data is generated by a standard process (such as a Poisson or a negative binomial 
in the case of accident frequencies, and a multinomial logit in the case of accident 
severities). Our presentation of Markov switching models is similar to that of Markov 



switching autoregressive models in econometrics 



McCuUoch and Tsayl . 



1994, 



Tsayl . 



2002 1 . 



All statistical models that we consider here, either for accident frequencies or for 
severity outcomes, are parametric and can be fully specified by a likelihood function 
f{Y\@,Ai), which is the conditional probability distribution of the vector of all 
observations Y, given the vector of all parameters of model Ai. If accident events 
are assumed to be independent, the likelihood function is 

T Nt 

/(Yie,A^)=nn^(^*."i®'-^)- (3.1) 



t=l n=l 



Here, Yt^n is the n observation during time period t, and P(F(,„|0,7W) is the prob- 
ability (likelihood) of Yt^n- The vector of observations Y = includes all obser- 
vations n = 1,2,. ..,Nt over all time periods t = 1,2, ...,T. Number Nt is the total 
number of observations during time period t, and T is the total number of time pe- 
riods. In the case of accident frequencies, observation 1^ „ is the number of accidents 
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observed on the n^^ roadway segment during time period t (note that Nt is the num- 
ber of roadway segments in this case). In the case of accident severity, observation 
Yf^n is the observed outcome of the n*^ accident occurred during time period t (note 
that Nf is the number of accidents in this case). Vector is the vector of all un- 
known model parameters to be estimated from accident data Y. We will specify the 
parameter vector separately for each statistical model presented below. Finally, 
model Ai = {M, ^t,n} includes the model's name M (e.g. M = "negative binomial" 
or M = "multinomial logit") and the vector ^t,n of all characteristic attributes (i.e. 
values of all explanatory variables in the model) that are associated with the n*^ 
observation during time period t. 



3.1 Standard count data models of accident frequencies 



The most popular count data models use d for predicting accident frequencies are 
Poisson and negative binomial (NB) models [Washington et al.l . l2003l |. These models 
are usually estimated by the maximum likelihood estimation (MLE) method, which 
is based on the maximization of the model likelihood function /(Y|0, A^) over the 
values of the model estimable parameters 0. 

Let the number of accidents observed on the n^^ roadway segment during time 
period t be At^n- Thus, our observations are Yt^n = A,™, where n = l,2,...,Nt and 
t = 1,2,...,T. Here Nt is the number of roadway segments observed during time 
period t, and T is the total number of time periods. The likelihood function for the 
Poisson m odel of accident frequenc ies is specified by equation (13. ip and the following 



equations 



Washington et al. 



2003j : 



nA,n\i3) 



Ay- 

\t,n = exp(/3'Xt,„), t 



exp(-At,„) 



1 2 T 



n 



l,2,...,iV,. 



(3.2) 

(3.3) 
(3.4) 
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Here, Xt^n is the Poisson accident rate for the n^^ roadway segment, this rate is equal 
to the average (mean) accident frequency on this segment over the time period t. 
The variance of a Poisson-distributed accident frequency is the same as its average 
and is equal to Xt^n- Parameter vector (3 consists of unknown model parameters to 
be estimated. Prime means transpose, so f3' is the transpose of (3. In the Poisson 
model the vector of all model parameters is @ = (3. Vector „ includes character- 
istic variables for the n^^ roadway segment during time period t. For example, Xj „ 
may include segment length, curve characteristics, grades, and pavement properties. 
Henceforth, the first component of vector Xj „ is chosen to be unity, and, therefore, 
the first component of vector f3 is the intercept. 

The likelihood function for the negative binomial (NB) model of accident frequen- 



cies i s specified by equation (13. ip and the following equations 



2003j : 



Washington et al. 



P{Yt,n\@,M) = P{At,n\e,M) = AfB{At,^\f3,a), (3.5) 
USiA..Ma, ^ Ii#f±iM f^V'° f^V-' . ,3.6) 



T{l/a)At^J \l + aXt,nJ \l + aXt,n, 
V = exp(/3%,„), t = l,2,...,T, n = l,2,...,Nt. (3.7) 

Here, r( ) is the standard gamma function. The over- dispersion parameter a > 
is an unknown model parameter to be estimated together with vector f3. Thus, the 
vector of all estimable parameters is = [/3', a]'. The average accident rate is equal 
to Xt^n, which is the same as in the case of the Poisson model. The variance of the 
accident rate is At,„(l + aXt^n), which is higher than in the case of the Poisson model 
(if a > 0). The negative binomial model reduces to the Poisson model in the limit 
a ^ 0. 

In addition to the Poisson and negative binomial models, we also consider the stan- 
dard zero-inflated Poisson (ZIP) and zero-inflated negative binomial (ZINB) models. 
These models account for a possibility of existence of two separate data-generating 
states: a normal count state and a zero-accident state. The normal state is unsafe, and 
accidents can occur in it. The zero- accident state is perfectly safe with no accidents 
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occurring in it 111 Zero-inflated models are usually used when there is a preponderance 
of zeros in the data. In the case of accident frequency data with many zeros in it, 
the probability of At^n accidents occurring on the n^^ roadway segment at time pe- 
riod t can be well modeled by a ZIP process or, if the data are over-dispersed, by a 



ZINB process. The likelihood functions of t 
equation (13.11) and the following equations 



le ZIP and ZINB mode^ 



Washington et al 



2003) : 



s are specified by 



P(V|0,A^) = P{At^n\@,M) 

= gt,„J(A,„) + (l-9t,n)P(A,n|/3) for ZIP, (3.8) 

P{Yt,n\@,M) = P{At,n\@,M) 

= qt,nIiAt,n) + il-qt,nWB{At,n\f3,a) for ZINB, (3.9) 



where 



J(A,n) = { 1 if A,n = 0, and if A,n > } , (3.10) 
1 

1 -\- e ^ 

Here we use two different specifications for the probability qt^n that the ra**^ road- 
way segment is in the zero-accident state during time period t. Scalar X^ n is the 
accident rate that is defined by equation (13. 4p . Probability distribution 2{At^n) is 
the probability mass function that reflects the fact that accidents never happen in 
the zero-accident state. The right-hand-side of equation (13. 8p is a mixture of the 
zero-accident distribution I{At^n) and the Poisson distribution V{At,n\/3) given by 
equation (13.31) . The right-hand-side of equation (13. 9p is a mixture of I{At^n) and the 
negative binomial distribution AfB{At^n\f3,a) given by equation (13.61) . Scalar r and 
vector 7 are estimable model parameters. We call "ZIP-r" and "ZINB-r" the models 
specified by equations (I3:8|) - (l3ll]) . We call "ZIP-7" and "ZINB-7" the models spec- 
ified by equations (13. 80 - 03. 101) and (13.120 . The vector of all estimable parameters is 



^Note that roadway segments are not required to stay in a particular state all the time and can 
move from normal count state to zero-accident state and vice versa. 
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e = [13', t]' for the ZIP-r model, = [/3',«,r]' for the ZINB-r model, = [/3',y]' 
for the ZIP-7 model, and = [/3',«,7']' for the ZINB-7 model. It is important to 
note that „ depends on the estimable model parameters and gives the probability 
of being in the zero-accident state, but qt^n is not an estimable parameter by itself. 



3.2 Standard multinomial logit model of accident severities 

The severity outcome of an accident is determined by the injury level sustained 
by the most severely injured individual (if any) involved into the accident. Thus, 
accident severities are a discrete outcome data. Most common statistical model used 
for predicting severity outcomes are the multinomial logit model and the ordered 
probit model. However, there are two potentia l problems with applying ordered 



Savolainen and Mannerind . 



2007 |. 



probability models to accident severity outcomes 
The first problem is due to under- reporting of non-injury accidents because they are 
less likely to be reported to authorities. This under-reporting can result in biased and 
inconsistent model coefficient estimates in an ordered probability model. In contrast. 



the coefficient estimates of an unordered mu 



for the intercept terms Washington et al. 



tinomial logit model are consistent except 



2003( 1 . The second problem is related to 



undesirable restrictio ns that ordered probabil ity models place on influences of the 



explanatory variables ^Washington et al. 



20031]. As a result, in this study we consider 



only multinomial logit models for accident severity. 

Let there be I discrete outcomes observed for accident severity (for example, 
J = 3 and these outcomes are fatality, injury and property damage only). Also 
let us introduce accident severity outcome dummies 5^^j^ that are equal to unity if 
the i^^ severity outcome is observed in the n^^ accident that occurs during time 
period t, and to zero otherwise. Then, our individual observations are the severity 
outcome dummies, Yt^n = {^t*!}, where i = 1,2,...,/. Note that n = 1,2,. ..,Nt and 
t = 1,2, ...,T, where Nf is the number of accidents observed during time period t, 
and T is the total number of time periods. The vector of all observations Y = {5^*^} 
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includes all outcomes observed in all accidents that occur during all time periods. The 
likelihood function for the multinomial logit (ML) model of accident severity outcorn es 
is specified by equation (13 .ip and the following equations 



Washington et al 



2003j : 



i=l 



i=l 

i 



1,2,..., I. 



(3.13) 
(3.14) 



Ej=iexp(/3jXt,„)' 

Parameter vectors /3j consist of unknown model parameters to be estimated, and 
f3 = {/3j}, where i = 1,2, I. Vector Xj „ contains all characteristic variables for the 
n^^ accident that occurs during time period t. For example, „ may include weather 
and environment conditions, vehicle and driver characteristics, roadway and pavement 
properties. We set the first component of X^ „ to unity, and, therefore, the first 
components of vectors /3j (z = 1, 2, /) are the intercepts. In addition, without loss of 
generality, we set all /3-parameters for the last severity outcome to zero, f3j = 0. This 
can be done without loss of generality because Xf ,^ are assumed to be independent 
of the outcome i, and, therefore, the numerator and denominator in equation (13. 14p 
can be multiplied by the an arbitrary common factor Washington et al.l . |2003 |. 



3.3 Markov switching process 

Let there be N roadway segments (or, more generally, roadway entities or/and 
geographical areas) that we observe during successive 

Markov switching models, which will be introduced below, assume that there is an 

unobserved (latent) state variable St^n that determines the state of roadway safety 

for the n^^ roadway segment (or roadway entity, or geographical area) during time 

period t. We assume that the state variable St^n can take on only two values: St^n = 

corresponds to the first state, and Sf^„ = 1 corresponds to the second state. The choice 

^In a more general case, we can observe a variable number of roadway segments over successive time 
periods. Here, for simplicity of the presentation, we do not consider this general case. However, our 
analysis is straightforward to extend to it. 



21 



of labels "0" and "1" for the two states is arbitrary and is a matter of convenience. 
We further assume that, for each roadway segment n the state variable St^n follows 
a stationary two-state Markov chain process in timelfl The Markov property means 
that the probability distribution of St+i^n depe nds only on the value St^n at time t, 



Breimanl . 



1969| . The stationary two- 



but not on the previous history St-2,', 
state Markov chain process {st,n} can be specified by time- independent transition 
probabilities as 

P{St+l,n = l\St,n = 0) = Po^^i, P{St+l,n = 0\St,n = 1) = p'Clo^ (3.15) 

where n = 1,2,...,N. In this equation, for example, P{st+i,n = l|sf,n = 0) is the 
conditional probability of St+i^n = 1 at time t + 1, given that St^n = at time t. 
Note that P(st+i,n = 0\st,n = 0) = = 1 - p'^q^.i and P{st+i,n = l|st,„ = 1) = 
p^^i = 1 — Pi'!lo- Transition probabilities Po"ii and p["io unknown parameters 
to be estimated from accident data (n = 1,2,...,A^). The stationary unconditional 
probabilities of states St^n = and St^n = 1 arc 

pin) ^ + p^'^p) for state Si,„ = 1. 

It is noteworthy that the case when (for each roadway segment n) the states St,n 
are independent and identically distributed in time t is a special case of the Markov 
chain process. Indeed, this case corresponds to history-independent probabihties 
of states "0" and "1", therefore, Po'!lo = Pi^o and Po*!!! = Thus, we have 

Pq^q = Pi"io — Po'^ ^^'^ Po'^i = Pi*^! = P\'\ where the last equalities in these two 
formulas follow from equations (13.161) . 



■^Stationarity of {sf_,i} is in the statistical sense [Breimanl . Il969l | 



^These can be found from the foUowing stationarity conditions: Pg"^ = [1 —pl^Xi\P<^'^ +?'i"!of^i"^ 
=p('^^pW + [l_p(n)^]p(») -{n)^.in)^^ [Breimanl . [l969t . 
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3.4 Markov switching count data models of annual accident frequencies 

When considering annual accident frequency data below, we will use and estimate 
two-state Markov switching Poisson (MSP) and two-state Markov switching negative 
binomial (MSNB) models that are proposed as follow. Similar to zero-inflated mod- 
els, these annual-accident-frequency Markov switching models assume that one of the 
two states of roadway safety is a zero-accident state, in which accidents never happen. 
The other state is assumed to be an unsafe state with possibly non-zero accidents oc- 
curring. MSP and MSNB models respectively assume Poisson and negative binomial 
(NB) data-generating processes in the unsafe state. Without loss of generality, below 
we take St^n = to be the zero-accident state and St^n = 1 to be the unsafe state. 

As in the case of the standard count data models of accident frequencies (see 
Section 13.11) . in this section, a single observation is the number of accidents At^n 
that occur on the ra**^ roadway segment during time period t. There are T time 
periods, each is equal to a year, and the periods are t = 1,2, ...,T. For simphcity o " 
presentation, we assume that the number of roadway segments is constant over time[f , 
Nt = N = const, and, therefore, the segments are n = 1,2, ...,N. The vector of all 
observations Y = {l^f,„} = {A,n} includes all accident frequencies At^n (t = 1,2, ...,T 
and n = 1,2, N). For each roadway segment n, the state St^n can change every 
year. The likelihood functions of the two-state Markov switching Poisson (MSP) and 
two-state Markov switching negative binomial (MSNB) models of annual accident 
frequencies A^^n are specified by equation (13.11) with A^^^ = A^, and by the following 
equations: 

f T(Atn) if = 

[ViAM itst,n = l 

^The analysis is easily extended to the case when we observe a variable number of roadway segments 
Nt 7^ const during time periods t, see also footnote [2] on page[20l In this case it would be convenient 
to count all segments asn=l,2,...,A^ and to count the time periods a,s t = T'y'^\ T^^^ + 1, T^"'', 

where the -n}^ segment is assumed to be observed during interval T^"''^ < t < rj"-* of successive time 
periods. 
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for the MSP model of annual accident frequencies, and 

f lUt n) if Stn = 

for the MSNB model of annual accident frequencies. Here zero-accident probability 
distribution X(y4j „), given by equation (13.101) . reflects the fact that accidents never 
happen in the zero-accident state St^n = 0. Probability distributions V{At^n\(3) and 
7Vi5(At^„|/3, a) are the standard Poisson and negative binomial probability mass func- 
tions, see equations (13. 3p and (13. 6p respectively. Vector f3 is the vector of estimable 
model parameters and a is the negative binomial over-dispersion parameter. To en- 
sure that a is non-negative, during model estimation we consider its logarithm instead 
of it. For each roadway segment n the state variable St^n follows a stationary two-state 
Markov chain process as described in Section 13.31 

Because the state variables „ are unobservable, the vector of all estimable pa- 
rameters must include all states (st,„), in addition to all model parameters (/?-s, 
a-s) and all transition probabilities (Po"^i, Pi^^q)- Thus, 

= [/3',a,p(l,...,pa,pS'-lo,-,pS'llsr, (3.19) 

where vector S = [{si^i, st,i), {si^n, st,n)]' contains all state values St^n and 
has length TxN. Of course, in the case of the MSP model, over-dispersion parameter 
a does not enter equation (13.191) . 

Note that, if Po^i < Pillo, then, according to equations (I3.16p . we have pg"'' > p^^\ 
and, on average, for the n^^ roadway segment state St^n = occurs more frequently 
than state St^n = 1- On the other hand, if p^f^i > Pi^lo, then state St^n = 1 occurs 
more frequently for the n*^ segment. 

In addition, note that here the choice of a year as the length of the time periods 
t = 1, 2, T is arbitrary. For example, one can consider quarterly (or other) periods 
instead. 

Finally, it is important to understand that although the MSP and MSNB models 
given by Equations (I3.17P and (13.181) assume state St^n = to be perfectly safe and 
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zero-accident, this state can be (and probably should be) viewed as an approximation 
for nearly safe states, in which accidents rarely occur. |^ 

3.5 Markov switching count data models of weekly accident frequencies 

When considering weekly accident frequency data below, we will use and estimate 
two-state Markov switching Poisson (MSP) and two-state Markov switching negative 
binomial (MSNB) models that are proposed as follows. In each of the two states 
(■St,n = and = 1) these weekly-accident-frequency models assume a standard 
Poisson data-generating process defined by equation (13. 3p . or a standard negative 
binomial process defined by equation (13.61) . Thus, both states are assumed to be 
unsafe for these models. We observe the number of accidents At^n that occur on 
the n^^ roadway segment during time period t, which is a week in this case. Let 
there be T weekly time periods in total. Let us again assume that the number 
of roadway segments is constant over time, Nt = N = const (see footnote [5] on 
page[22]). Thus, in equation (13. ip the vector of all observations is Y = {l^,n} = {^t,n}, 
where t = 1,2, ...,T and n = 1,2,. ..,N. In addition, for weekly-accident-frequency 
Markov switching models, we assume that all roadway segments always have the 
same state, and, therefore, the state variable St^n = St depends on time period t only. 
This is because, here, state St is intended to capture common unobserved factors 
influencing roadway safety on all segments. Correspondingly, all roadway segments 
switch between the states with the same transition probabilities Po^li = Po^i and 

(n) 

With this, the likelihood functions for the two-state Markov switching Poisson 

(MSP) and two-state Markov switching negative binomial (MSNB) models of weekly 

^Nearly safe states have average accident rates At.„ <C 1 [see Equations (|3.4p and (|3.7p ]. In this case, 
the perfectly safe, zero-accident state, which has At,„ = 0, serves as a good approximation for these 
nearly safe states. 
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accident frequencies At^n are specified by equation (13.11) with A^^ = A^, and by the 
following equations: 

P{Yt^n\®,M) = P{At.n\@,M) = < 

[ P(A,n|/3(i)) if Si 
for the MSP model of weekly accident frequencies, and 

T^n^ ir^ ^,^ T^,A\r^KA\ 1 ■^^( An|/3(o), «(0)) 

\ Ari3(A,n|/3(l),«(l)) 

for the MSNB model of weekly accident frequencies. Here, t = 1,2,...,T and n = 
1,2,...,N. Probability distributions V{. . .) and AfB{. . .) are the standard Poisson 
and negative binomial probability mass functions, see equations (13. 3p and (13. 6p re- 
spectively. Parameter vectors /3(q) and /3(i), and negative binomial over- dispersion 
parameters a(o) > and q;(i) > are the unknown estimable model parameters in 
the two states St = and St = I. To ensure that a(o) and q;(i) are non- negative, 
their logarithms are considered during model estimation. Because, we choose the 
first component of Xt^„ to be equal to unity, the first components of /3(o) and f3(^i^ 
are the intercepts in the two states. Note that the state variable St follows a station- 
ary two-state Markov chain process with transition probabilities po^i and pi-,o as 
described in Section [3731 

Because the state variables St are unobservable, the vector of all estimable param- 
eters must include all states (st), in addition to all model parameters (/5-s, a-s) 
and all transition probabilities {po^i, Pi^o)- Thus, 

0= [/3(o),tt{o),/3(i),a(i),Po-.i,Pi->o,S']'- (3.22) 

where vector S = [si,...,st]' has length T and contains all state values. In the 
case of the MSP model, over- dispersion parameters a(o) and are absent from 
equation (I3.22p . 




1 



(3.20) 



if St = 



(3.21) 
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Without loss of generality, we assume that (on average) state St = occurs more 
or equally frequently than state St = 1. Therefore, po > pi, and from Equations (I3.16P 
we obtain restrictiorQ 

Po^i < Pi^o- (3.23) 

In this case, we can refer to states St = and Si = 1 as "more frequent" and "less 
frequent" states respectively. 

Note that here the choice of a week as the length of the time periods t = 1, 2, T 
is arbitrary. For example, one can consider daily (or other) periods instead. 



3.6 Markov switching multinomial logit models of accident severities 

When considering accident severity data below, we will use and estimate two- 
state Markov switching multinomial logit (MSML) model that is proposed as follows. 
In each of the two states (0 and 1), this model assumes standard multinomial logit 
(ML) data-generating process that is defined by equation fl3.14p and described in 
Section 13.21 We observe severity outcome dummies that are equal to unity if the 
i^^ severity outcome is observed in the n^^ accident that occurs during time period t, 
and to zero otherwise. We consider weekly time periods, t = 1,2,...,T, where T is 
the total number of periods observed. Then, the vector of all observations Y = {5^^} 
includes all outcomes observed in all accidents that occur during all time periods, 
i = 1,2,..., I, n = l,2,...,Nt and t = 1,2, ...,T. Here I is the total number of 
possible severity outcomes, and Nt is the number of accidents observed during weekly 
time period t. For MSML models of accident severities, we again assume that all 
roadway segments (where accidents happen) always have the same state of roadway 
safety, and, therefore, the state variable St^n = St depends on time period t only (in 
this case, state St captures common unobserved factors that infiuence safety on all 

^Restriction (I3.23|) is introduced for the purpose of avoiding the problem of switching of state labels, 
1. This problem would otherwise arise because of the symmetry of the likelihood functions 
given by equations (13. ip . (|3.20p and p.2ip under the label switching. 
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segments). Correspondingly, all roadway segments switch between the states with 
the same transition probabilities pl^lli = Po^i and p^iIIq = Pi^o- 

The likelihood function for the two-state Markov switching multinomial logit 
(MSML) model of accident severity outcomes is specified by equation (13.11) and the 
following equations: 

P(V|0,A^) = l[[P{i\@,M)f 



n5« 



1 = 1 



n [X/:(z|/3(o))]''" if3t = 



i=l 



n [^£(^1/3(1))]'"" iist = l 

i=l 



(3.24) 



where n = l,2,...,Nt and t = 1,2, ...,T. Probability distributions A^£(z|/3(o)) and 
A4C{i\l3(^i-^) are the standard multinomial logit probability mass functions in the two 
states, see equation (13.141) . Here /3(q) = {/3(o),j} and = where i = 

1, 2, /. Parameter vectors /3(q) j and /3(i),j are unknown estimable model parameters 
in states and 1 respectively. Since we choose the first component of X^^^ to be equal 
to unity, the first components of vectors /3^q-^ ^ and /3(^) ^ are the intercepts. Similar to 
the case of the standard (single-state) ML model presented in Section 13.21 here, we 
set all /5-parameters for the last severity outcome to zero, /3(o),/ = /^(i),/ = 0- 

The vector of all estimable parameters includes all states (st), in addition to 
all model parameters (/9-s) and all transition probabilities (po^i) Pi^o)- Thus, 

= [/3;o),/3;,),Po^i,Pi^o,Sr. (3.25) 

where vector S = [si, s^]' has length T and contains all state values. 

In analogy with the assumption made in the previous section, here, without loss of 
generality, we assume that (on average) state = occurs more or equally frequently 
than state St = 1. Therefore, po > Pi, and from equations (I3.16P we again obtain 
restriction 

Po^i < Pi^o- (3.26) 
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In this case, we can refer to states St — and = 1 as "more frequent" and "less 
frequent" states respectively. 

Note that here the choice of a week as the length of the time periods t = 1, 2, T 
is arbitrary. For example, one can consider daily (or other) periods instead. 
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CHAPTER 4. MODEL ESTIMATION AND COMPARISON 



This chapter presents the basics of Bayesian estimation of standard models and 
Markov switching models of accident frequencies and severities. We also discuss com- 
parison of different models by using Bayesian approach, and an evaluation of model 
fit performance. 



4.1 Bayesian inference and Bayes formula 

Statistical estimation of Markov switching models is complicated by unobservabil- 
ity of the state variables St^n (or St)^ As a result, the traditional maximum likelihood 
estimation (MLE) procedure is of very limited use for Markov switching models. 
Instead, a Bayesian inference approach is used. Given a model Ai with likelihood 
function f(Y\@,M.), the Bayes formula is 

/(0|Y ^), /(Y.e|A^) ^ /(Y|e.xMe|A^) 

Here /(0|Y,A^) is the posterior probability distribution of model parameters 
conditional on the observed data Y and model Ai. Function f{Y,@\A4) is the 
joint probability distribution of Y and given model A4. Function f{Y\Ai) is the 
marginal likelihood function - the probability distribution of data Y given model Ai. 
Function 7r(0|A^) is the prior probability distribution of parameters that reflects prior 
knowledge about 0. The intuition behind equation (14. ip is straightforward: given 
model Ai, the posterior distribution accounts for both the observations Y and our 



^For example, in the case of Markov switching models of weekly accident frequencies, we will have 
260 time periods (T — 260 weeks of available data). In this case, there are 2^^° possible combinations 
for value of vector S — [si, st]'- 



30 



prior knowledge of 0. We use the harmonic mean formula to calculate the marginal 



likelihood /(Y|A^) of data Y [see iKass and Rafteryl . Il995| as 



f{Y\M)-' = /(Y|A<)-'|,r(e|7M)<ie = /(Y|M)-'|^||j|M<i0 



/(Y|0,A1) 



where E{. . . |Y) is the posterior expectation (which is calculated by using the posterior 
distribution). 

In our study (and in most practical studies), the direct application of equa- 
tion (14.11) is not feasible because the parameter vector contains too many com- 
ponents, making integration over in equation (14. ip extremely difficult (see foot- 
note[T]on pagel^H]). However, the posterior distribution /(0|Y, A^) in equation (14. ip 
is known up to its normalization constant, namely /(0|Y,A^) oc f{Y,@\Ai) = 
f(Y\@,M)7r{@\M). As a result, we use Markov Chain Monte Carlo (MCMC) sim- 
ulations, which provide a convenient and practical computational methodology for 
sampling from a probability distribution known up to a constant (the posterior dis- 
tribution in our case). Given a large enough posterior sample of parameter vector 
0, any posterior expectation and variance can be found and Bayesian inference can 
be readily applied. In the next chapter we describe our choice of prior distribution 
7r(0|A^) and the MCMC simulations in detail. The prior distribution is chosen to be 
wide and essentially noninformative. For the MCMC simulations, we wrote a special 
numerical code in the MATLAB programming language and tested it (for details see 
the next chapter). 

In the end of this section, let us make a short noteworthy digression. In Bayesian 
statistics model observations and model parameters are treated on an equal footing. 
Therefore, for Markov switching models, one can treat the vector of all state values 
S as latent model parameters, or as latent (hidden) observations. We treat S as 
model parameters. As a result, in our approach, the transition probabilities Pi"^>o 
and Poj'>i do not enter the likelihood function /(Y|0, A^), which is a function of 
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S and model coefficients (/3-s, as, 7-s, r) only (refer to the likelihood functions 
presented in the previous chapter). In this case, the Markov switching property is 
treated as a prior information, and the prior distribution, given in the next chapter, 
reflects this property (in other words, we a priori specify that the state variable St^n 
follows a Markov process in time). If we treated state values S as latent observations, 
then the vector of all observation would include both Y and S. In this case, the 
likelihood function would depend on the transition probabilities and would become 
/(Y,S|e\S,7W) = f{Y\&,M)f{S\&\S,M), where 0\S means all components of 
except S. In any case, for the purpose of model comparison discussed below, the 
marginal likelihood should always be defined as /(Y|A^) [not as /(Y, S|A^)] because 
Y is the only data that is truly observed. 



4.2 Comparison of statistical models 

For comparison of different models we use the following Bayesian approach. Let 

there be two models Aii and Ai2 with parameter vectors ©i and ©2 respectively. 

Assuming that we have equal preferences of these models, their prior probabilities are 

7r(A^i) = 7r(A^2) = 1/2. In this case, the ratio of the models' posterior probabilities, 

P{A4i\Y) and P(A^2|Y), is equal to th e Bayes factor. The late r is defined as the 

ratio of the models' marginal likelihoods [Kass and Rafteryl . Il995i] . Thus, we have 

P{M2\Y) ^ f{M2,Y)/f{Y) ^ f{Y\M2)7r{M2) ^ fjYlM^) . . 

PiMi\Y) /(A^i,Y)//(Y) /(Y|^i)7r(A^i) /(Y|A<i)' ^ ' ' 

where f{Aii,Y) and f{Ai2,Y) are the joint distributions of the models and the 
data, /(Y) is the unconditional distribution of the data, and the marginal likelihoods 
/(Y|A^i) and /(YIAI2) are given by equation (14. 2p . If the ratio in equation (14.31) is 
larger than one, then model Ai2 is favored, if the ratio is less than one, then model A4i 
is favored. An advantage of the use of Bayes factors is that it has an inherent penalty 
for including too many parameters in the model and guards against overfittingy 

^There arc other frequently used model comparison criteria, for example, the deviance information 
criterion, DIG = 2E[D{@)\Y] - D{E[&\Y]), where deviance £1(0) = -2 ln[/(Y|0, 7W)] [Roberl, 
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4.3 Model performance evaluation 



To evaluate the performance of m odel {.M,©} in fitting; the observed data Y 



we carry ou t a goodness-o: 



Wood . 



2002 



Press et al. 



2007 |. 



-fit test 



Maher and Summersgilll . 



1996 



Cowanl . 



1998 . 



In the case of accident frequency models, quantity X 



is! 



X 



T Nt 

EE 

t=l n=l 









&,M) 



(4.4) 



where E(Yt^n\®, -M) and var(Yt^n\&, -M.) are the expectations and variances of the 
observations Yt^n- In accident frequency studies, the observations are the frequencies, 
Yt^n = At-n on roadway segment n during time period t. For example, from equa- 
tions (13. 6p . (13. 7p and (13.211) for the MSNB model of weekly accident frequencies we 
find the following formulas for the (unconditional of state) expectations and variances: 



E{Yt^^\@,M) = p^Xfl+p^Xfl and mr(Fi,„|0, ^) = poAi;:^ (1 + «(o)Ai;:^) +PiAi;^(l + 



"(1)^*3 +PoPi(aJ11 - \^t,nf, where X'fl = exp(/3(o)Xt,„) and X\ll = exp(/3(i)Xt,„) are 
the mean accident rates in the states St = and St = 1 respectively. For the MSNB 
model of annual accident frequencies one needs to set Xf^l = in these formulas be- 
cause state = is the zero-accident state in this case. The appropriate formulas 
for Poisson models can be obtained by setting the over- dispersion parameters (as) 
to zero. 

In the limit of asymptotically normal distribution of large accident frequencies, 
has the chi-square distribution with degrees of fr eedom eg ual to the number of 
observations minus the number of model parameters Wood . |2002| |. Because weekly 
(and even annual) accident frequencies are typically small, in this study, we do not 
rely on the assumption of their asymptotic normality. Instead, we carry out Monte 



,(1) 



l200l| . Models with smaller DIG are favored to models with larger DIG. However, DIG is theoreti- 
cally based on the assumption o f asymptotic multivariate normality of the posterior distribution, in 
which case DIG reduces to AIG jSpiegelhalter et aL . 2002| . As a result, we prefer to rely on a math- 
ematically rigorous and formal Bayes factor approach to model selection, as given by equation (|4.3p . 
■^Note that for a standard Poisson distribution, the variances are equal to the means, 
var{Yt_n\®,J^) ~ E{Yt^n\®,-^), and equation (|4.4p reduces to the Pearson's x^- 
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Carlo simulations to find the distribution of [Cowanl . |1998| . This is done by 
generating a large number of artificial data sets under the hypothesis that the model 
{Ai, 0} is true, computing and recording the value for each data set, and then 
using these values to find the distribution of x^- This distribution is then used to 
find the goodness-of-fit p-value, equal to the probability that exceeds the observed 
value of (the later is calculated by using the observed data Y) Q 

In the case of accident severity models, we use the Pearson's x^, defined as 

T Nt I 



^ 2-2-2-. P{z\@,M) ' ^ ' 

t=l n=l 1=1 \ I ' / 

where the accident severity outcome dummies d^^l are equal to unity if the i*^ severity 
outcome is observed in the n**^ accident that occurs during time period t, and to zero 
otherwise. According to equation fl3.24p . the theoretical unconditional probability of 



the i^^ outcome is P{i\&,M) = poMC{i\f3^Q-j) + piMC{i\l3 



^Note that for this Monte Carlo simulations approach, specification of quantity is actually very 
flexible. For example, one can potentially use [Yt^n~E{Yt^n\®, ■M)]'^ /var{Yt^n\®, under the sum 
in equation ()4.4|) for the goodness-of-fit test. However, in this case would not become chi-square 
distributed even in the asymptotic limit of large accident frequencies. 
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CHAPTER 5. MARKOV CHAIN MONTE CARLO SIMULATION 

METHODS 



We use MCMC simulations for Bayesian inference and model estimation. This chap- 
ter presents MCMC simulation methods in detail. First, we describe a hybrid Gibbs 
sampler and the Metropolis-Hasting algorithm. Next, we explain a general Markov 
switching model representation that we use for all Markov switching models of acci- 
dent frequencies and severities. After that we describe our choice of prior probability 
distribution. Then we give detailed step-by-step algorithm used for our MCMC sim- 
ulations. Finally, in the end of this chapter, we briefly overview several important 
computational issues and optimizations that allow us to make Bayesian-MCMC es- 
timation reliable, efficient and numerically accurate. For brevity, in this chapter we 
omit model specification notation in all equations. For example, in this chapter 
we write the posterior distribution /(0|Y, A^) simply as /(0|Y), and etc. 



5.1 Hybrid Gibbs sampler and Metropolis-Hasting algorithm 

As we have mentioned in the previous chapter, because the posterior distribution, 
given by the Bayes formula fl4.1l) . is extremely difficult to find exactly, but is relatively 
easy to find with accuracy up to its normalization constant, we use Markov Chain 
Monte Carlo (MCMC) simulations. They provide a feasible statistical methodology 
for sampling from any probability distribution known up to a constant, the posterior 
distribution in our case. 

To obtain draws of the parameters vector from a posterior distribution /(0|Y), 
we use the hybrid Gibbs sampler, which is an MCMC simulation algorithm that 
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involves both Gibbs and Metropolis -Hasting; sampling McCuUoch and Tsayl . Il994 . 



Tsay 



2002 



SAS Institute Inc. 



2006| . Assume that is composed of K components: 
= [6[,62, ...,6'j^y , where 6^ can be scalars or vectors, k = 1,2,...,K. Then, the 
hybrid Gibbs sampler works as follows: 

1. Choose an arbitrary initial value of the parameter vector, = 0'^'^) , such that 
/(0W|Y) > [i.e. /(0(°)|Y) oc /(Y, 0(o)) = /(Y|0W)7r(0(o)) > 0]. 

2. For each g = 1,2,3,... , parameter vector 0*^^^ is generated component-by- 
component from 0(9-1) by the following procedure: 

(a) First, draw 6[^^ from the conditional posterior probability distribution 
f{0^f^\Y, O'f ^\ 9^^^^^). If this distribution is exactly known in a closed 
analytical form, then we draw 0^^^ directly from it. This is Gibbs sampling. 
If the conditional posterior distribution is known up to an unknown nor- 
malization constant, then we draw 6^f^ by using the Metropolis-Hasting 
(M-H) algorithm described below. This is M-H sampling. 

(b) Second, for all k = 2, 3, — 1, draw 0^^^ from the conditional posterior 
distribution f{ei^^\Y, e[^\ ^^^2^, e^j^^l\ e^f^^) by using either Gibbs 
sampling (if the distribution is known exactly) or M-H sampling (if the 
distribution is known up to a constant). 

(c) Finally, draw 6^^^ from the conditional posterior probability distribution 
/(0^''|Y, 9^^\ 0^^Li) by using either Gibbs or M-H sampling. 

3. The resulting Markov chain {@^^^} converges to the true posterior distribution 
/(0|Y) as g oo. 

Note that all conditional posterior distributions are proportional to the joint distri- 
bution /(Y, 0) = /(Y|0)7r(0). For example, we have 

fro o a a a \ f(y,0i,---,0k-i,0k,0k+i,---,0K) 

J{tfk\^ , ffi, ffk-i, ffk+i, tiK) - a a a a \ 

J[Y ,tfl, Uk-l, ffc+l, •••5 f^K) 

cx /(Y, 01, dk, 6>fc+i, Ok) = /(Y, 0). (5.1) 
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By using the hybrid Gibbs sampler algorithm described above, we obtain a Markov 
chain {0(^)}, where g = 1, 2, ... , Gbi, Gu + 1, . . . , G. We discard the first Gu "burn- 
in" draws because they can depend on the initial choice Of the remaining 
G — Gbi draws, we typically store every third or every tenth draw in the computer 
memory. We use these draws for Bayesian inference. We typically choose G ranging 
from 3 X 10^ to 3 x 10^, and Gm = G/10. In our study, a single MCMC simulation run 
takes from one day to couple weeks on a single computer CPU. We usually use eight 
different choices of the initial parameter vector Thus, we obtain eight Markov 



chains of 0, and use them for the Brooks-Gelman-Ru 
of our MCMC simulations 



Brooks and Gelman 



jin diagnostic of convergence 



19981], for details see Section 15.51 



below. We also check convergence by monitoring the likelihood f(Y\@^^^) and the 
joint distribution f(Y,@^^^). 

We use the Metropolis-Hasting (M-H) algorithm to sample from conditional pos- 
terior distributions known up to their normalization constantsll Specifically, our goal 
here is to draw e[^^ from f{ek\Y,d[^\...,o[^l^,o[^~l\...,e^^~^^) distribution that is 
not known exactly, so we cannot use the Gibbs sampling. The M-H algorithm works 
as follows: 

• Choose a jumping probability distribution J{0k\0k) of Ok- It must stay the 
same for all draws g = Gbi + 1, G, and we discuss its choice below. 



Draw a candidate 0k from J{6k\0 



Calculate ratio 



(9-1) ^ 



^ ""1(9-1) IV aid) a(9) fl(9-l) /3(9-1)n TtO \q(9-^)\' ^ 



k 



Set 



(g) _ I Ok with probability min(p, 1), 

^k — \ yp-^) 

6): otherwise. 



^In general, the M-H algorithm allows to make draws from any probability distribution known up 
to a constant. The algorithm converges as the number of draws goes to infinity. 
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Note that the unknown normahzation constant of /(. . .) cancels out in equation (15.21) . 
Also, if the jumping distribution is symmetric J{0k\0k) = J{0k\6k), then the ratio 
J{0l^ k) I J {b k\6^^ ^'*) becomes equal to unity and Metropolis-Hasting algorithm 
reduces to Metropolis algorithm. The averaged acceptance rate of candidate values 
in equation ( 15. 3p is recommended to range from 15 to 50%. In this study, during the 
first Gbi burn-in draws we make adjustments to the jumping probability distribution 
j(0k\6k) in order to achieve a 30% averaged acceptance rate during the Metropolis- 
Hasting sampling (carried out during the remaining G — Gu draws used for Bayesian 
inference). The specifics about the choice of the jumping distribution and of its 
adjustments are given below in Sections 15.41 - 15. 5[ 

5.2 A general representation of Markov switching models 

All Markov switching models for accident frequencies and severities, specified in 
Sections 13.41 - [3l6l can be represented in a general, unified way. This representation 
allows us to estimate all models by using the same mathematical notations, compu- 
tational methods and, most important, the same numerical code. In this section, 
first, we introduce a convenient general representation of Markov switching models 
considered in this study. Second, we show how Markov switching models for accident 
frequencies and severities, specified in Sections 13.41 - 13.61 are described by using this 
general representation. 

For the general, unified representation of Markov switching between the roadway 
safety states over time, we would like to make the state variable to be dependent 
on time only. For this purpose, we introduce an auxiliary time index t, so that the 
state variable si depends only on i. For example, in the case of annual frequencies of 
accidents occurring on roadway segments over T annual time periods (this case is 
considered in Section [3^ . the auxiliary time is defined as t = t + (n — 1)T, where the 
real time is t = 1,2,...,T and the roadway segment number is n = 1,2,...,A^. The 
auxiliary time index runs from one to N x T, that is t = 1,2, ...,NT. As another 
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t: 1 2 3 4 5 6 7 8 9 10 11 



^ ^' l^o->r l^i->o ^ ^' l^o->r 



(r=2) „(r=2) 
l->0 



Figure 5.1. Auxiliary time indexing of observations for a general Markov 
switching process representation. 



example, consider the case of weekly accident frequencies observed over T weekly 

time periods (refer to Section [3.51) . In this case the auxiliary time simply coincides 

with the real time, t = t. 

A general scenario of Markov switching between the roadway safety states over 

auxiliary time t is schematically demonstrated in Figure 15. 1[ The auxiliary time 

index runs from one to T, that is i = 1,2, ...,T. During an auxiliary time period i 

the system is in state (which can be or 1). As the auxiliary time index increases 

from t to t + 1, the state of roadway safety switches from to s^^i. We assume that 

for all t ^ 71 (for all t that do not belong set 71) this switching is Markovian, that 

is the probability distribution of Sf_^i depends on the value of (see Section 1X51) . 

We assume that for those values of t that belong to the set 71, the switching is 

independent of the previous state, that is for t G 71 the probability distribution of 

s,^, is independent of and of the earlier states]] The values t G T_ are shown 

by white dots in Figure (15. ip . the values t ^ are shown by black dots, and the 

Markov switching transitions are shown by concave arrows. In a general case, the 

transition probabilities for Markov switching si s^^^, where t ^ 71, do not need to 

be necessarily constant and can depend on the auxiliary time index t. As a result, we 

assume that there are R auxiliary time intervals T{r) < t < T{r + 1), r = 1, 2, i?, 

^Independent switching can be view as a special case of Markovian switching, see the discussion that 
follows equation (|3.16p 
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such that the transition probabihties are constant inside each time interval and can 
differ from one interval to another. Here the set T contains, in an increasing order, 
all left boundaries of the time intervals, the first element of T is equal to 1, and 
the last element of T is equal to T + 1. Note that the size of set T (i.e. the 
number of elements in it) is equal to i? + 1. Thus, to repeat, for each value of index 
r = 1,2, the transition probabilities Po'2,i and p^'^g constant inside the r^^ 
interval T(r) < t < T{r + 1). In Figure (15.11) the intervals of constant transition 
probabilities are shown by curly brackets beneath the dots. 

In the real time t all data observations (accident frequencies or severity outcomes) 
are counted by using the real time index, that is the vector of all observations is 

Y = {yi,n}, where t = 1,2,...,T and n = 1,2, Nt. When we change to the 
auxiliary time, all observations are counted by using the auxiliary time index, that is 

Y = {yjfj}, where t = 1, 2, T and n = 1, 2, iV^. Here A^^ and iV^ are the number 
of observations during real and auxiliary time periods t and t respectively. There is 
always a unique correspondence between the indexing pairs (t,n) and (t, n). Using 
the auxihary time indexing, the likehhood function /(Y|0), given by equation (13.11) . 
becomes 



f Ni 



i=l n=l 



t "="i [ /(>f,ril3(i)) if si = 1 



Nr 



{i: .31=0} n=l 



Nr 



{i: 31=1} n=l 



(5.4) 



where /(^t,ji|/3(o)) and /(^t,ji|/3(i)) are the model likelihoods of single observations 
.j^ in roadway safety states = and = 1 respectively. Set {i: = 0} is defined 
as all values of t such that 1 < t < T and = 0, and set {t : = 1} is defined 
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analogously. Vectors /3(-o) and are the model parameters vectors in the states 
and 1, these vectors are specified by the model type as follows: 



/3 



/3(-^-) for Poisson or multinomial logit, 

[f3'u),a{s)]' for negative binomial, 

(5.5) 



Hs)]' or [/3(^), a(^), r(,)]' for ZIP-r or ZlNB-r, 
, Wis)n[s)]' or 7(.)]' for ZIP-7 or ZINB-7 models, 

where s = 0, 1 are the state values. Scalar r and vector 7 are estimable zero-infiated 
model parameters, and a is the over-dispersion parameter, as defined in Section 13. 1[ 
By defining the auxiliary time t and sets 7^ and T, we specify the general unified 
representation of the Markov switching models introduced in Chapter [31 as follows: 

• For Markov switching models of annual accident frequencies, introduced in Sec- 
tion [3]4l we have 

i=t + {n~l)T, f = NxT, n = l, iV^ = 1, (5.6) 
71 = {nT, where n = 1, A^}, (5.7) 
T = {l + {r-l)T, (l + NT)}, r = l,...,N, R = N, (5.8) 
n = \i/T] and t = i - {n - 1)T, (5.9) 

where t = 1, 2, T and n = 1,2, ...,N are the real time index and the roadway 
segment number respectively, and [a;] is the "ceil" function that returns the 
smallest integer not less than x. Here T is the number of annual time periods, 
and is the number of roadway segments observed during each period. The 
change of indexing to auxiliary time t, given by equation (15.61) . is demonstrated 
in Figure 15.11 for the case when T = 5 (in Section 16.11 we will consider a five- 
year accident frequency data). Separate roadway segments n = 1,2,..,N have 
different transition probabilities for their states of roadway safety [refer to equa- 
tion (13.151) ]. Therefore, in Equation (15.81) the time interval number r coincides 
with the roadway segment number n, that is r = n and R = N. Equation (15. 7p 
follows from the fact that states Sf are independent for different roadway seg- 
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ments n = 1,2,.. .,N. Equation (15.91) gives the conversion from the auxihary 
time indexing back to the real time indexing. 

The observations are annual accident frequencies „ (refer to Sections 13.11 
and 13.41) . Therefore, we have Y^f^ = Yi^ = Yt-n = At^n, where t and n are 
calculated from i by using equations (15. 9p . Thus, according to equations (I3.17P 
and (I3.18p . the likelihood functions of a single observation l^^^ = Y^^^ = At^n in 
the states and 1 are 

/(n^ii^(o)) = /(n,ii^{o)) = AAt,n). 
f{yiSm) = /(nii^ci)) = ^(Ani^d)) 

for the MSP model of annual accident frequencies, and 

for the MSNB model of annual accident frequencies. Here n = 1, while t and n 
are calculated from t by using equations (15. 9p . Keep in mind that /3(]^) is given 
by equation (15.51) . 

For Markov switching models of weekly accident frequencies, introduced in Sec- 
tion [3i5l we have 

i=t, f = T, n = n, Ni = N, (5.12) 
r_ = {0}, T = {1, (T + 1)}, r = l, R=l, (5.13) 

where t and n are the real time index and roadway segment number respectively, 
T is the number of weekly time periods, and is the number of roadway 
segments observed (it is the same for all periods). Here the auxiliary time t 
coincides with the real time t. The transition probabilities are constant over all 
time periods t = t and are the same for all roadway segments n = 1,2, ...,N. 
Thus, R = 1, set T consists of just two values, and set 71 is empty. 

The observations are weekly accident frequencies At^n (refer to Section 13.50 . 
Therefore, we have = Yt^n = ^t.n, where we use i = t and h = n. Thus, 
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according to equations fl3.20p and fl3.2ip . the likelihood functions of a single 
observation F^^^ = A^^n in the states and 1 are 

/(%lAo)) = ^(A,„|3(0)), /(n^^l^d)) = ^(A,n|3(l)) (5.14) 

for the MSP model of weekly accident frequencies, and 

/(n,nl^(0)) = Ari3(A,n|^(0)), /(^nl^d)) = Ari3( A,n |^(1) ) (5.15) 

for the MSNB model of weekly accident frequencies. Here t = t and n = n. 
Note that /3(q) and /3(-]^) are given by equation (15.51) . 



For Markov switching models of accident severities, introduced in Section | 
we again consider weekly time periods and, therefore, have formulas very similar 
to equations (l532|) -( l5:T3ll . 

t = t, f = T, h = n, Ni = Nt, (5.16) 
T_ = {0}, T = {1, (T + 1)}, r = l, R=l. (5.17) 

Here, the auxiliary time t again coincides with the real time t, scalar T is 
the total number of weekly time periods, and Nt is the number of accidents 
occurring during time period t. 

The observations are accident severity outcome dummies 5^*^ (refer to Sec- 
tion [XS]). Thus, we have Y^f^ = Ytn = {^t*l}, where i = 1,2,..., I and we 
use t = t and n = n. According to equation fl3.24p . the likelihood functions of 
a single observation If ^ in the states and 1 are 



i=l 

f{Yi,Sa)) = n [-^^(^13(1))]''" , (5.18) 



where t = t and n = h. Note that /3(o) and /S^j^^ are given by equation (15. 5p . 
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In the remaining sections of this chapter we use the above general representation 
of Markov switching models. For convenience and brevity of the presentation, we 
drop tildes (~) from all our notations. In other words, we use t, T, n, Nt and f3 
instead of t, T, n, iV^ and /3. We also call "auxiliary time" just "time". Thus, it 
is good to keep in mind that, in the rest of this chapter, time index/period/interval 
means auxiliary time index/period/interval. 



5.3 Choice of the prior probability distribution 

A full specification of Bayesian methodology and model estimation requires a 
specification of the prior probability distribution. In this section we describe how we 
choose the prior distribution 7r(0) of the vector of all parameters to be estimated. 
In our study, for the general representation given in the previous section, vector 
includes all unobservable state variables (st), model parameters (/3(q), and 
transition probabilities for every r^^ time interval (Po^i, Pi^l^o, = 1,2, ...,R). Thus, 

= [/3;o),/3(i),rf-li,---,pSi,p£o,---,pSo,ST- (5.19) 

Here, vectors /3(q-) and /3(-^-) are the model parameter vectors for states s = and 
s = 1, which are defined in equation fl5.5l) . Vector S = [si, S2, st]' contains all 
state values and has length T, which is the total number of time periods. 

The pri or distribution is suppose d to refiect our prior knowledge of the model 



SAS Institute Inc. 



2006| . We choose the prior distributions of /3(qv, /3 



(1)' 



parameters 

PqI^i and p^{\q (r = 1, 2, i?) to be nearly fiat and essentially non-informativeif| The 
prior distribution of the state vector S must refiect the Markov switching property of 
the state variable St. The overall prior distribution of the vector of all parameters 
is chosen to be the product of the prior distributions of all its components [refer to 
equation fl5.19p ]. Thus, our choice of the prior is as follows: 



■^equation ()4.1|) shows that for nearly flat prior distributions, when 7r(0|A^) is approximately con- 
stant around the peak of the likelihood function, the posterior distribution only weakly depends on 
the exact choice of the prior. We have verified this result during our test MCMC runs. 
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• Prior probability distribution of model parameters vectors is the product 
of prior distributions for the vector components in states s = and s = 1, 

1 ^(s) 

7r(/3(o),/3(i)) = nn^(A^)''^)' (^-20) 

s=0 k=l 

where P{s),k is the k^^ component of vector /3(g), and K(^s) is the length of vector 
/3(-g-) (i.e. the number of model parameters in the state s is equal to -^'(s), where 
s = 0, 1). For free parameters (3(s),k (which are free to be estimated), the priors of 
P{s),k are chosen to be normal distributions: 7!'{P{s),k) = J^{l3{s),k\fJ'k, Sfc). [Keep 
in mind that for NB models ln(a) is estimated instead of the over- dispersion 
parameter a, and, thus, the prior distribution of a is log-normal.] Parameters 
that enter the prior distributions are called hyper-parameters. For these, the 
means /i^ are chosen to be equal to the maximum likelihood estimation (MLE) 
values of Pk for the corresponding standard single-state models (Poisson, NB, 
ZIP, ZINB and multinomial logit models in this study). The variances are 
chosen to be ten times larger than the maximum between the MLE values of 
jSk squared and the MLE variances of jSk for the corresponding standard models 
(thus, variances are chosen to be relatively large in order to have wide prior 
distributions of P(s),k)- 

All /3-parameters can be either free (which are free to be estimated) or restricted 
(which are not free to be estimated, but instead are set to some predetermined 
values). We choose normally-distributed priors only for free parameters. In this 
study, if a parameter is not free, then there are only three other possibilities: the 
non-free parameter is restricted to be equal to either zero, or — oo, or a free pa- 
rameter. Thus, in all these three cases we have prior knowledge about the value 
of the restricted parameter. For simplicity of presentation, in equation fl5.20p 
and below we do not explicitly show which /3-parameters are free and which are 
restricted, and for presentation purposes only we portray all /5-parameters as 
being free. However, it is important to remember that during numerical MCMC 
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simulations we do not draw restricted parameters 
the appropriate values that they are restricted to|^ 



Dut, instead, we set them to 



For weekly accident frequency and severity models, introduced in Sections 13.51 
and ESI the joint prior distribution for all transition probabilities {pq}^i,Pi}^q}, 
where r = 1, 2, i? (note that i? = 1 in case of basic weekly models), is 

vr({rili,Pr-lo}) oc fl7r{ptlMptl,)I{pi;l, < pZo)- (5.21) 



r=l 



Here T^iPol^i) = Betalp^l^ilvojUo) and T^ip^ 



(r) 



JO 



Beta{p^{}^Q\vi,i'i) are cho- 



sen to be standard beta distributions. Function I{pq}^^ < p^^q) is defined as 
equal to unity if restriction p^}^i < p^{}^o is satisfied and to zero otherwise [re- 
fer to equation fl3.23l) ]. For annual accident frequency models, introduced in 
Sections 13.41 the prior distribution for transition probabilities is given by equa- 
tion fl5.2ip with functions I{pq}^i < Pi^l^o) dropped out because there are no any 
restrictions for transition probabilities in this case [note that equation (15.211) be- 
comes an equality in this case]. Thus, in the case of annual accident frequency 
models, functions I{p^}^i < Pi'l.o) should be left out from all formulas in the 
rest of this chapter. The hyper-parameters in equation (15.211) are chosen to be 
■^0 = '^o = "^1 = '^i = 1) in which case the beta distributions become the uniform 
distribution between zero and one. Similar to parameters I3[s),ki we draw only 
free transition probability parameters Pq\,i and p^i\,Q- All restricted transition 
probabilities are not drawn, but are set to the values that they are restricted 
to. 



non-free parameter that is restricted to a free parameter is set immediately after the free param- 
eter is drawn during the hybrid Gibbs sampler simulations. This is because these two parameters 
(the restricted "child" parameter and its "parent" free parameter) must always be the same. For 
example, if we have three beta-parameters /3i, j32 and /^a, and if (3^ is restricted to then (3^ is set 
to the new value of /3i immediately after this new value is drawn. 



46 



The prior distribution for the state vector S = [si, S2, st]' is equal to the 
hkehhood function of S given the transitional probabilities {pq}^i,p^{}^q}, 

/(siui,p!io})=n^i) n n 

(t:l<t<T,\ Jt:l<t<T,\ 

( t: l<t<T, \ 

R 

=n n pist+i\st) 

r = l r t: T(r)<t<T{r + l),\ 
\ t<T, t^T_ j 

= fl [p^olir^^ [1 - ptlr^^^ i&r^^^^ [i - p^^lof^^^ ■ (5.22) 

r=l 

Here, index r = 1, 2, R counts time intervals T(r) < t < T(r + 1) of constant 
transition probabilities Po^i and p^{}^q (see Section r5.2p . Number m[^2,j is the 
total number of Markov switching state transitions from St = i to St+i = j 
inside time interval T{r) < t < T{r + 1) [here i,j = {0,1} and history- 
independent transitions at t G 71 are not counted]. In equation (15.221) we 
disregard probability distribution P{si) and distributions P{st+i) for t G 71 . 
This is because their contribution is negligible when T is large and the number 
of elements in set 71 is small relative to the value of T, which is true in this 
studylf] 

It is important to note that formula (15.221) for distribution /(S|{po^^,p['2,o}) 
follows from the Markov switching property of the state variable St [compare 
to equation (13.151) ]. In other words, we a priori specify that the state variable 
St follows a Markov process in time, with transition probabilities {pq1,i, p^{1,q} , 
and this specification must be and is reflected in the prior distribution given by 
equation fl5.22p . 



5 



Alternatively, we can assume that P(si = 0) = P(si = 1) = 1/2 and P{st+i = 0) = P{st+i — 
1) = 1/2 for all t e 71. A more sophisticated alternative (not considered here) would be to treat 
these probabilities as free estimable parameters of the model. 
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• Finally, the prior probability distribution 7r(0) of parameter vector 0, which 
is given by equation (15.191) . is the product of the priors of all 0's components, 
given by equations (15.201) - (I5.22p . 

7r(e) = 7r(S,Ul„rilo},/3(o),/3(i)) 

= /(S|U'ii,pSlo})^(U':li,pS'io})^(/3(o),/3(i)) 

( t: l<t<T, \ 

L t^T- S 
B. 

X JJi3eta(p[,'2,i|i;o, uq) Beta{pf IqIvi, ui)I{pi[l^ < p^Iq) 

r=l 
1 

xllll^fiP(s),kW,^k). (5.23) 

s=0 k=l 



5.4 MCMC simulations: step-by-step algorithm 

In our research, for Bayesian inference on the parameter vector 0, given by equa- 
tion (I5.19p . we apply the hybrid Gibbs sampler and make draws of the components 
of vector from their conditional posterior distributions (refer to Section [STT]) . All 
conditional posterior distributions are proportional to the joint distribution /(Y, 0) 
[see equation (15.11) ]. The joint distribution is /(Y, 0) = /(Y|0)7r(0), where the 
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likelihood /(Y|0) is given by equation (15.41) and the prior 71(0) is given by equa- 
tion (I5.23p . Thus, for the joint distribution we have 



/(Y,e) = /(Y|0)7r(0 

Nt 

OC 



Nt 



{t:st=l}n=l 



X 



n n/(^*."i/5(o)) 

{t: st=0} n=l 

n pi't+M 

( t: l<t<T, \ 

R 



r=l 



X 



K, 



(0) 



(1) 



JjA/'(/3(l),fc|yUfc,S^ 



(5.24) 



As a result, the conditional posterior distributions of all components of vector 0, 
which are proportional to the joint distribution, are as follows: 

• The conditional posterior distribution of the k^^ component of vector f3(Q) is 



/(Ao),fc|Y, 0\/5(o),fc) = ^ — Y— ^ /(Y' ®) 



OC 



Nt 



n n/(^*."i/5(o)) 

{f : st=0} n=l 



/(Y,0\/3(o),,) 

X 7V(/3(o),fc|/U/c, Sfc) 



Aft 



{t: st=0} n=l 



1 



X 



OC 



Nt 



n n/(vi/3(o)) 

{<: st=0} n=l 



X g-[/3(0),fc-Mfc]V2Sfc_ 



-[/3(0),fc— Mfel^/SSfc 



(5.25) 



where 0\/?(o),a; means all components of except /?(o),fe, and we keep only those 
multipliers that depend on /?(o),a:- In equation (15.251) the conditional posterior 
distribution of /5{o),fc is known up to an unknown normalization constant. There- 
fore, we draw free parameters /5(o),fc by using the Metropolis-Hasting algorithm 
described in Section [5TT1 Note that k = 1,2, ...,K(^q), where i^'(o) is the number 
of model's /3-parameters in state 0. 
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The conditional posterior distribution of the k component of vector /3(i), 
is derived similarly to the conditional posterior distribution of P(o),k in equa- 
tion fl5:25|) . 



/(/3(i),fc|Y,0\/5(i),,)cx/(Y,0) 



oc 



Nt 



{t:st=l}n=l 



(5.26) 



Free parameters P(i),k, where k = 1,2, K(^i), are also drawn by using the 
Metropolis-Hasting algorithm. 

(r) 

The conditional posterior distribution of the transition probability Po-*i 



/(p£i|Y,0UZi 



(r) 



(r) 



f(Y, 0\p(,l, 



oc/(Y,0) 



oc Y\_ Pist+i\st) 

( t: l<t<T, \ 



,(0 



f=l 
X 



OC 



i [Vo)L [Uo) 



oc Seta(p[|'2,i|m[|'2,i + t^o, "^oZo + '^o)/(Po-ii < Pi 



(r) 



(5.27) 



where r(. . .) is the Gamma function, 0\pq'2,i means all components of ex- 



cept Po'2»i, and we keep only those multipliers that depend on PqI^i- We use 
formula (15.221) to obtain the fourth line in equation (I5.27p . and number mf^j is 
the total number of Markov switching state transitions from st = i to St+i = j 
inside time interval T(r) < t < T{r + 1) [not counting history-independent 
transitions at t G 71]. In equation (I5.27P the conditional posterior distribution 
of Pq\i is a known truncated beta distribution. Therefore, we draw Pq\,i di- 
rectly from this distribution by using Gibbs sampling described in Section 15.11 



(r) 
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and by the rejection sampling ("accept-reject") algorithm described in the next 
section. Note that index r = 1,2, ...,R, where R is the total number of time 
intervals of constant transition probabilities. 

• The conditional posterior distribution of the transition probability p^{}^o is given 
by equation fl5.27p with states and 1 interchanged everjrwhere, except in func- 
tion liPoii < Pilo)^ 

/(pSlo|Y,e\p(lo)oc/(Y,0) 

cx i3eta(p[io|mSlo + Vi,m^;l, + Pi)I{ptli < pSIo)- (5-28) 

We also draw Pil_,Q directly from its conditional posterior distribution by using 
Gibbs sampling. 

• To speed up MCMC convergence for posterior draws of vector S = [si, S2, st]', 
we draw subsections St^r = [st, St+i, st+r-i]' of S at a time, instead of draw- 
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ing vector S component-by-component Tsayl . |2002| . The conditional posterior 
distribution of is 



/(Y,ew) 



oc 



X 



{i: Sf=0} "=1 



X 



n n/(nni/3(i)) 

{t: S£=l} »i=l 



t: l<t<T, 

iir_ 



OC 



n n/(nj/3(o)) 



X 



t<t<t + T--l 

n n 



n n/(nj/3(i)) 



t<t<t + T — 1 



r=l 



t+r-1 

n 



t: T(r)<t<T(r + l), t<T, 
t-l<f<t + T-l, f^T_ 



n=l 



n=l 



X 



n ipiiii'<"!ii -p£j*"ipnoi*"ii -pr-ioi*-' 



r=l 



t+r-1 



n 



Afr 



?i=l n=l 

X n uiii*-' [1 - fiir- b'S'-ioi"- [1 - j>n„i*-' , (5.29) 

{v. [r{r),r(r+l))n [t-l,t+r-l]^{0}} 

where 0\Sf means all components of except for and we keep only those 
multipliers that depend on St^,- = [sj, St+i, St+r-^ ■ Number is the total 
number of Markov switching state transitions from S£ = i to = j inside the 
intersection of time intervals T(r) < i < T{r + 1) and t — 1 <i<t + T — l [here 
= {0)1} ^'^d history- independent transitions at t G 71 are not counted]. 
Number rhi^j is zero for all i,j = {0,1} if intervals T{r) < i < T{r + 1) 
and t — l<£<t + r — Ido not intersect, resulting in the final expression for 
the product over r on the last line in equation fl5.29p . Vector St^,- has length 
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r and can assume 2^ possible values. By choosing r small enough, we can 
compute the right-hand-side of equation (15.291) for each of these values and find 
the normalization constant of /(St^T-|Y, 0\Sf^^). This allows us to make Gibbs 
sampling of St^r- Out typical choice of r is from 5 to 14. 

All components of parameter vector are given by equation (I5.19p . and all con- 
ditional posterior distributions are given by equations fl5.25l) - fl5.29l) . We generate 
draws of ©'•^•' from 'j^y using the hybrid Gibbs sampler explained in Section ISTTl 

as follows (for brevity, we drop g indexing below): 

(a) We draw vector /3(-q-) component-by-component by using the Metropolis-Hasting 
(M-H) algorithm (note that we draw only those components that are free pa- 
rameters). For each (free) component /3(o),fe of /3(q) we use a normal jumping 
distribution 

AMPio),k)=Af0iomo),k,al^,,) = L=e-I'^(°)''=-^<«)-'=lVH),. (5.30) 

O"(o),fcv27r 

Standard deviations <J(o),k are adjusted during the burn-in sampling (i.e. during 
g = 1,2, ...,Ghi) to have approximately 30% averaged acceptance rate in equa- 
tion (15.31) . The adjustment algorithm is explained in the next section. We also 
tried Cauchy jumping distribution 

J0{o),k\(^io),k) = Cauchy0(o)\P{o),k,crio),k) 

l/(7r(T(o),fc) 



1 + 



(Ao),fc - /5(o),fc)Mo),fc 



2 ' 



(5.31) 



and obtained similar results. As already explained in Section [5731 we draw P{o),k 
from its conditional posterior distribution, given by equation (I5.25p . only if it 
is a free parameter. We do not draw /?(o),fc in the following three cases. First, 
/3(o),fc is restricted to zero (which is the case if it is found to be statistically 
insignificant). Second, /3(o),fc is restricted to —oo [which is the case if state 
is the zero-accident state, and, therefore, the intercept in state is —oo, 
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see equations (13.41) . (13. 7p . (13.171) and (13.181) ]. The third case is when /?(o),fe is 
restricted to another, free /3-coefficient. 

(b) We use Metropohs-Hasting algorithm and draw all components of (3(^i^ (that are 
free parameters) from their conditional posterior distributions, given in equa- 
tion (I5.26p . in exactly the same way as we draw the components of /3(o)- 

(c) By using Gibbs sampling, for all r = 1,2, ...,R time intervals we draw transi- 
tion probabilities Pq1,i, first, and Pil^Q, second, from their conditional posterior 



distributions given in equations (15.271) and (15.281) 



(d) Finally, we draw subsections St^r = [st, St+i, St+r-i]' of the state vector S = 
[si, S2, St]' ■ We use Gibbs sampling and draw subsections St^r one after 
another from their conditional posterior distributions given by equation (15.291) . 



5.5 Computational issues and optimization 



A special numerical code was written in the MATLAB programming language for 
the MCMC simulations used in the present research study. Our code was written from 
scratch, and no standard MCMC computer scripts and procedures were used. This 
programming approach provided us with great fiexibility and control during model 
estimation. Our code uses the general representation introduced Section 15. 2[ and as 
a result, the code is applicable to estimation of all accident frequency and severity 
models considered here. 

^We do not make draws of Pg'^,! and p^iXq from their conditional posterior distributions if these 
parameters are not free, but are restricted to other transition probabihties. For example, in the 
next chapter we will consider a model of weekly accident frequencies in which we will assume that 
different seasons have different transition probabilities, but the transition probabilities for the same 
seasons in different years are restricted to be the same. In this case, only transition probabilities for 
time intervals that are inside the first year of data are free and are drawn. 
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Below, in this section, we briefly discuss several numerical issues, tips and opti- 
mizations that turned out to be important for numerically accurate, reliable and fast 
MCMC runs during model estimation process. 

• We tested our MCMC code on artiflcial accident data sets. The test procedure 
included a generation of artiflcial data with a known probabilistic model (e.g. 
a MSNB model or a MSML model). Then these data were used to estimate the 
underlying model by means of our MCMC simulation code. With this procedure 
we found that the probabilistic models, used to generate the artiflcial data, were 
reproduced successfully with our estimation code. 

• In order to avoid numerical zero and numerical infinity, during MCMC simula- 
tions we always use and calculate the logarithms of all probability distributions 
instead of the distributions themselves (for example, we work with log-likelihood 
functions instead of likelihood functions). 

• Standard deviations cr(o),k of the normal and Cauchy jump distributions, given 
by equations (15.301) and (15.311) . are adjusted during the burn-in sampling {g = 
1,2, ... ,Ghi) to have approximately 30% averaged acceptance rate in equa- 
tion (15.31) . For each k = 1, 2, ...-ft'(o) that corresponds to a free model parameter 
P(o),k, drawn by the Metropolis-Hasting (M-H) algorithm, the adjustment is 
done as follows. We calculate the mean candidate acceptance rate in equa- 
tion (15.31) . averaged over the last 50 consecutive M-H draws. If this mean rate 
is above/below the 30% target rate, we multiply/divide the standard deviation 
"^(o),*,' by factor 1.25. Then we calculate the mean acceptance rate, averaged over 
the next 50 M-H draws, and again adjust <7(o)^k by multiplying or dividing it by 
1.25, and so on. During the burn-in sampling we collect and save all standard 
deviations used for the M-H draws and the corresponding mean acceptance rates 
(averaged over groups of 50 consecutive draws). After all Gbi burn-in draws are 
made, we fit a decreasing exponential function to the dependence of the mean 
acceptance rates on the cr(o),fc values [for this flt we use the acceptance rate data 
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collected over the last {2/3)Gbi burn-in draws]. Finally, we use this exponential 
function to obtain the best guess about the value of o"(o),a,. that will result in the 
30% target averaged acceptance rate. This value of cr(o),fc stays constant for all 
further draws g = + 1, ...,(?, which are used for Bayesian inference. As a 
result of the adjustments descried above, in our MCMC simulations the actual 
mean acceptance rates (averaged over draws g = Ghi + 1, always turned 

out to be within 1% of the target 30% rate. 

The Gibbs sampling draws from the truncated betas distributions in equa- 
tions (I5.27P and (15.281) are don e by the reiection samp ling technique, also known 



as the accept-reject algorithm Hormann et al.l . |2004| . This algorithm works as 



follows. Let us assume that we need to make draws of x from a probabil- 
ity density function /(x), which is not easily available. Then, we construct 
an envelope function F{x) such that, first, F{x) > f{x) is satisfied for all 
X, and, second, x can be easily drawn from the probability density function 
F{x) j J F{x) dx. To obtain correct draws from f{x), we repeatedly, first, gen- 
erate draws Xg from F(x) j J F{x) dx, and, second, accept Xg with probability 
f{xg)/F{xg) [here g = 1,2,3,...]. For the algorithm to be efficient, the enve- 
lope function F{x) should be sufficiently close to f{x) [so that the acceptance 
probability f{xg)/F{xg) is not very small]. Because the logarithm of a trun- 
cated beta distribution is concave, we construct and use a piece-exponential 



env elope function 
see Hormann et al 



i .e. th e logarithm of the envelope function is piece-linear). 



2004 1 



The Gibbs sampling of subsections St^r = [st, St+i, St+r-i]' from the condi- 
tional posterior distribution given in equation fl5.29p can be optimized as follows. 
First, for each value of time t = t, t + 1, t + r — 1 we calculate and save in the 
computer memory the values of products Yl f0^i,n\(^{o)) Yl f(Xi,n\(^{i))^ 

n=l n=l 

refer to equation (15.291) . Then, we use these values to compute the probabilities 
of all 2^ possible combination values of the subsection vector St^r without need 
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to recalculate the likelihood functions /(^£,„|/3(o)) and each time. 

This optimization procedure considerably speeds up Gibbs samphng of St^r- 

There is an important issue that arises during Bayesian-MCMC estimation of 
Markov switching models, which is the "label switching problem". This prob- 
lem can be understood and solved as follows. Note that the likelihood func- 
tions for the MSP, MSNB and MSML models, given by equations fl3:20D . 
and (13.241) . are completely symmetric under the interchange "0" of the la- 

bels of the two states of roadway safety. This label interchange is just equivalent 
to renaming labels for the two states (using label names "1" and "0" as opposed 
to using label names "0" and "1" for the first and second states respectively). 
During a MCMC run the labels might interchange many times back and forth, 
in which case the MCMC chain would not converge. This is called the "label 
switching problem" . To avoid this problem, we impose a restriction po-^i ^ Pi^o 
on the Markov transition probabilities, see equations (13.231) and (I3.26p . This 
restriction breaks the symmetry of the likelihood function and the posterior dis- 
tribution under the interchange "0"^"1" of the label notations|j In practice, 
the restriction imposed on the transitional probabilities does not completely 
solve the label switching problem because few MCMC chains still happen to 
converge to the incorrect label setting (with the two labels interchanged as 
compared to the correct label setting). To deal with this problem, we moni- 
tor the posterior average of the logarithm of the joint probability distribution 
/(Y, 0) [note that monitoring the joint distribution is equivalent to monitoring 
the posterior distribution because the later is proportional to the former] . When 
a MCMC chain converges to an incorrect label setting, this average is consider- 
ably smaller (typically, by 10 to 50) than its value for the MCMC chains that 



Instead of the restriction imposed on the transitional probabilities, we also tried restrictions imposed 
on the intercept coefficients (the first components of (3-a). We found that the later works no better 
and no worse that the former for controlling the label switching problem. It is convenient to use 
the restriction on the transitional probabilities because there are more than two intercepts in the 
MSML models and because of its easier interpretation [the interpretation of restriction po->i ^ Pi—^o 
is that, on average, the state is more frequent than the state 1, refer to equation l|3.16p ]. 
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converge to the correct label setting. To distinguish label settings, we define the 
correct label setting as the one that provides the maximal value of the average 
of the posterior probability and, therefore, the maximal value of the average 
of the joint probability (since the posterior is proportional to the joint). If we 
had an unlimited computational time, then eventually all MCMC chains would 
converge to the correct label setting. Since our computational time is limited, 
we have to eliminate those few chains that did not converge to the correct label 
settings!^ 

The convergence of the MCMC chains used for Bayesian inference (these chains 
have converged to the correct label setting) is checked by the Br o oks- G elman- 



Rubin diagnostic Brooks and Gelmanl . Il998l . ISAS Institute Inc.l . |2006| |. This 



diagnostic works as follows. Let us consider only the continuous model param- 
eters (/5-s, Po-ii and PiI^q). Let the number of these continuous parameters 
be equal to and vector be composed of these parameters (in other words, 
d include all components of except the state values S). Let us have M 
chains with G posterior draws of in each chain, obtained as a result of the 
MCMC simulations^ Thus, we have draws 0^3, ^ where g = 1,2,...,G and 
m = 1,2,...,M. Then, the potential scale reduction factors (PSRFs) and the 



8 



This may introduce a model estimation bias. However, this bias is neghgible because the incorrect 
label setting corresponds to posterior (or joint) probability values that are much smaller than those 
for the correct label setting (typically the difference factors range from « e^^" to « e~^°). 
^We usually use eight MCMC chains of vector of all model parameters, which correspond to 
different choices of the initial parameter vector ©'"^ see Section [STTl Keep in mind that a certain 
number of first "burn-in" draws of are discarded in each chain because they can depend on the 
initial choices of 0'"'. 
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multivariate potential scale reduction factor (MPSRF) are given by the follow- 
ing equations: 

G M 

g=l m=l 

M 



m=l 

^ M r- G 

w = — Y -J—Y(e^3'^^-e^'^\e 



{g,m) ^(m)^^^(c/,m) _ g{m)y 



m=l 3=1 



PSRF. = ^Ip^^, where n = l,2, N, (5.32) 
diag(Vr)„ 



MPSRF = + ^^±^X^,,, (5.33) 

where Xmax = max {eigenvalues of matrix W^^B}. 

Here vectors 6^^''^\ o'^'^^ and 6 have their lengths equal to A^; matrices 5, W 
and V have their sizes equal to N xN; PSRF„, is the potential scale reduction 
factor for the n^^ component of and diag(\^)„ is the n^^ element of the diagonal 
of matrix V [diag(iy)„ is defined analogously]; Xmax is the maximal eigenvalue 
of the symmetric positive definite matrix W^^B. Matrices B and W represent 
variance/covariance of between the MCMC chains and inside the chains, re- 
spectively. For a well-converged MCMC simulation, the resulting PSRFs and 
MPSRF should be close to unity. Note that there exist an alternative definition 
of PSRF and MPSRF, in which the square roots are missing in Equations fl5.32p 
and ^Mj- 
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CHAPTER 6. FREQUENCY MODEL ESTIMATION RESULTS 



In this chapter we present model estimation results for accident frequencies. The 
chapter consists of two sections. In the first section, we consider annual accident 
frequencies and estimate Markov switching Poisson (MSP), Markov switching nega- 
tive binomial (MSNB), standard Poisson, standard negative binomial (NB), standard 
zero-inflated Poisson (ZIP) and standard zero-inflated negative binomial (ZINB) mod- 
els. In the second section, we consider weekly accident frequencies and estimate MSP, 
MSNB, standard Poisson and standard NB models. We compare the performance of 
the models in fitting the data. 

In the present study, for both annual and weekly accident frequency models, we use 
the data from 5769 accidents that were observed on 335 interstate highway segments 
in Indiana in 1995-1999. 



6.1 Model estimation results for annual frequency data 

In this section we use annual time periods, t = 1, 2, 3, 4, T = 5 in totalQ For each 
roadway segment n = 1,2, . . . , N = 335 the state can change every year. Three 
types of annual accident frequency models are estimated: 

1. We estimate standard (single-state) Poisson and negative binomial (NB) mod- 
els, specified by equations (13. 3p and (13. 6p . We estimate these models, first, 
by the maximum likelihood estimation (MLE) and, second, by the Bayesian 

^We also considered quarterly time periods and obtained qualitatively similar results (not reported 
here) . 
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inference approach and MCMC simulations^ As one expects, for our choice 
of a non-informative prior distribution, for both the Poisson and NB models, 
the estimated results obtained by MLE and by MCMC estimation techniques, 
turned out to be very similar. 

2. We estimate standard zero-inflated ZIP-r, ZIP-7, ZlNB-r and ZINB-7 models, 
specified by equations fl3.8l) - fl3.12l) . First, we estimate these models by the 
MLE (see footnote [2] on page EH]). Second, we estimate them by the Bayesian 
inference approach and MCMC simulations. As one expects, for our choice of 
a non-informative prior distribution, the Bayesian-MCMC estimation results 
again turned out to be similar to the MLE estimation results for the ZIP-r and 
ZINB-r models. 

3. We estimate the two-state Markov switching Poisson (MSP) and two-state 
Markov switching negative binomial (MSNB) models, given in equations f l3.17p 
and fl3.18p . by the Bayesian-MCMC methods. To choose the explanatory vari- 
ables for the final MSP and MSNB models reported here, first, we start with 
using the variables that enter the standard Poisson and NB models^ Then, we 
consecutively construct and use 60%, 85% and 95% Bayesian credible intervals 
for evaluation of the statistical significance of each /^-coefficient in the MSP 
and MSNB models. As a result, in the final MSP and MSNB models some 
_n. of , „e ..need .o .e J Fo. MSNB .ode... no ..nc.on . 



2 



The maximum likelihood estimation was done by using LIMDEP software package. To obtain 
optimal parsimonious standard models, estimated by MLE, we ch oose the explanatory variables 



and their dummies b y using the Akaike I nformation Criterion ( AIC) [Tsavl . I2002L [Washington et al 
E003]. For details see iMalvshkiTal |2006l |. 



■^This approach makes comparison of explanatory variable effects in different models straightforward. 
A formal Bayesian approach to model variable selection is based on evaluation of model's marginal 



likelihood and the Bayes factor (j4.3p . Unfortunately, because MCMC simulations are computation- 
ally expensive, evaluation of marginal likelihoods for a large number of trial models is not feasible 
in our study. 

^A /3-coefhcient is restricted to zero if it is statistically insignificant. A 1 — a credible interval is 
chosen in such a way that the posterior probabilities of being below and above it are both equal to 
a/2 (we use significance levels a = 40%, 15%, 5%). 
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imposed on the over-dispersion parameter a, which turns out to be statistically 
significant anyway. 

The estimation results for the standard Poisson and NB models of annual accident 
frequencies are given in Table 16.11 The estimation results for the zero-infiated and 
Markov switching Poisson models are given in Table 16.21 The estimation results 
for the zero-infiated and Markov switching negative binomial models are given in 
Table 16.31 In these tables, posterior (or MLE) estimates of all continuous model 
parameters, (3-s and a, are given together with their 95% confidence intervals (if 
MLE) or 95% credible intervals (if Bayesian-MCMC), refer to the superscript and 

Table 16.41 gives 

summary statistics of all roadway segment characteristic variables „ except the 
intercept. 

Because estimation results for Poisson models are very similar to estimation results 
for negative binomial models, let us focus on and discuss only the estimation results for 
negative binomial models. Our major findings, discussed below for negative binomial 
models, hold for Poisson models as well (unless otherwise stated). The findings are 
as follows. 

The estimation results show that two states of roadway safety exist, and that 

the two-state MSNB model is strongly favored by the empirical data, as compared 

to the standard ZINB-r and ZINB-7 models, which in turn are favored over the 

simple standard NB model. Indeed, from Tables 16.11 and 16.31 we see that the values 

of the logarithm of the marginal likelihood of the data for NB, ZINB-r, ZINB-7 and 

MSNB models are -2554.16, -2519.90, -2447.33 and -2184.21 respectively. Thus, 

the MSNB model provides considerable, 369.95, 335.69 and 263.12, improvements of 

the logarithm of the marginal likelihood as compared to the NB, ZINB-r and ZINB-7 

models respectively. As a result, from equation (14. 3p . we find that, given the accident 

^Note that MLE assumes asymptotic normality of the estimates, resuhing in confidence intervals 
being symmetric around the means (a 95% confidence interval is ±1.96 standard deviations around 
the mean). In contrast, Bayesian estimation does not require this assumption, and posterior distri- 
butions of parameters and Bayesian credible intervals are usually non-symmetric. 
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data, the posterior probability of the MSNB model is larger than the probabilities of 
the NB, ZINB-r and ZINB-7 models by e^^^-^^, e^^^-^^ and e^^s-ia respectivelyO Note 
that we use the harmonic mean formula, given in equation fl4.2p . to calculate the 
values and the 95% confidence intervals of the log-marginal-likelihoods reported in 
Tables I6.1H6.3[ The confidence intervals are found by 

We can also use a classical statistics approach for model comparison, based on the 
MLE. Referring to Tables 16.11 and I6.3[ the MLE gives the maximum log-likelihood 
values -2533.81, -2502.67 and -2426.54 for the NB, ZINB-r and ZINB-7 models 
respectively. The maximum log-likelihood value observed during our MCMC simu- 
lations for the MSNB model is equal to —2049.45. An imaginary MLE, at its con- 
vergence, would give MSNB log-likelihood value that would be even larger than this 
observed value. Therefore, if estimated by the MLE, the MSNB model would provide 
large, at least 484.36, 453.22 and 377.09, improvements in the maximum log-likelihood 
value over the NB, ZINB-r and ZINB-7 models. These improvements would come 
with no increase or a decrease in the number of free continuous model parameters 
(/?-s, a, r, 7-s) that enter the likelihood function. Both the Akaike Information Cri- 
terion (AIC) and the Bayesian Information Criterion (BIC) would strongly favor the 
MSNB models over the NB modelEI 



^In addition, we find DIG (deviance information criterion) values 5105.4, 5037.3, 4891.4, 4261.5 for 
the NB, ZINB-r, ZINB-7 and MSNB models respectively. We also find DIG values 5346.0, 5292.9, 
5063.8, 4317.0 for the Poisson, ZIP-r, ZIP-7 and MSP models respectively. This means that the 
MSNB (MSP) model is favored over the standard NB (Poisson) and ZINB (ZIP) models. However, 
we prefer to rely on the Bayes factor approach instead of the DIG (see footnote [2] on page[3T|). 
^During bootstrap simulations we repeatedly draw, with replacement, posterior values of to 
calculate the posterior expectation in equation ()4.2|) . In each of 10^ bootstrap draws that we make, 
the number of values drawn is 1/100 of the total number of all posterior values available from 
MCMC simulations. The bootstrap simulations show that equation (|4.2p gives sufficiently accurate 
answers, and that the expectation in this equation is not dominated too much by just few posterior 
values of 0, at which the likelihood function happens to be extremely small. 

^Minimization of AIC ~ 2K — 2LL and BIC = K\ n(N) — 2LL ensures an opt imal choice of 
explanatory variables in a model and avoids overfitting Tsav . 20021 Washington et al. . 2003], Here 
K is the number of free continuous model parameters that enter the likelihood function, N is the 
number of observations and LL is the log-likelihood. When N > 8, BIG favors fewer free parameters 
than AIC does. 



Table 6.1 

Estimation results for standard Poisson and negative binomial models of annual 



accident frequencies 



Variable 


Poisson 


NB 


by MLE^ 


by MCMC ^ 


by MLE-^ 


by MCMC^ 


Intercept (constant term) 


-15.71 


-15.71 


-20.0Zlti 


— 20.3_23 1 


Accident occurring on interstates 1-70 or 1-164 (dummy) 


-.689I:-«« 


-.6891:??° 


-•756::«;]« 


-•76o::gg 


Pavement quality index (PQI) average® 


— .023o 


-.oi84--!"i51 

— .0234 


-.0150-°°^f 

— .0235 


-.0149-°°^f 

— .0231 


Road segment lengtli (in miles) 




.0504-!!IP!^ 

.025(3 






Logarithm of road segment length (in miles) 


.924 ggg 


•925-i?? 

",871 




.ggo^dSI 

'-^'^".938 


Number of ramps on the viewing side 


-.0397-°"? 


.0396_ Qg^g 






Number of ramps on the viewing side per lane per mile 


.335 


.414i°2 

.335 


.407f?i 

' .312 


.4ioi}2 

".312 


Number of lanes on a roadway 






.513I}/; 




Median configuration is depressed (dummy) 


177.277 
' ' .0775 


1 7©. 278 
' ".0788 


^ ' .0558 


^^.0o4o 


Median barrier presence (dummy) 


'-'•'-'^—3.63 


— 3.06 

— 3.71 


-2.41-2" 
— 2.82 


-2.44-V^? 

— 3.02 


Interior shoulder presence (dummy) 


-1.09-,^?t 

— i.YO 


— i.81 






Width of the interior shoulder is less that 5 feet (dummy) 


•358:ii 


Qc:a-467 
.aoo 


.OJO 207 


.358:|;t 


Outside shoulder width (in feet) 




_ 0614- °^^° 


— Ofi'^2-'°2®^ 


.UD0.3_ gg^g 


Outside barrier absence (dummy) 


-244i:Ji 


-■244I:Jg 


- 251--"l 


nc:n-.0984 
■^"'^-.406 


Average annual daily traffic (AADT) 


-3.981^:15 
X 10-5 


— ■^ qq— 3.17 

X 10^5 


-4.83:;::^^ 

X 10"" 


—4 qq-3-84 
^■=''^-6.08 

X lO-"' 


Logarithm of average annual daily traffic 




2.05?il 


9 9c;2.61 


oq2.()1 
^•^°1.96 


Posted speed limit (in mph) 


0121 O^""' 


0121 -"204 

.UiZi „g37g 


■'-'-'-^>'. 000762 


Q-|4f;.0280 
■'J -"-^".00128 


Number of bridges per mile 


00860 


- 0262" 

.UZDZ_ g43g 


-026ii:°r7r 


m7n--"0652 
— •'-'^"-'-.0488 


Maximum of reciprocal values of horizontal curve radii (in 1 / mile) 


-1641:1°^ 


— .165_ 222 


-.i94i:^°f 


— IQfi— '^^2 
.i»u_ 282 


Maximum absolute value of change in grade of a vertical curve 


■0456:0iii 


■0456:«i| 






Number of vertical curves per roadway section 




-158::*^ 






Percentage of single unit trucks (daily average) 


1 401-*'' 


l-40^9!i 


i-ee^sio 


1.662gti 


Number of changes per vertical profile along a roadway segment 


nfiifiio^ 

.UDiD gj^4Q 


■0619;io°i 


.0597;i«l9 


•06i4:j;ii 



Table 6.1: (Continued) 



Variable 


Poisson 


NB 


by MLE'^ 


by MCMC'' 


by MLE<^ 


by MCMC<^ 


Ovcr-([ispc'rsiuii ])araiiic'lc'r fi in XJJ iiiu<[cls 










Mean accident rate {Xt,n for Poisson and NB), averaged over all values of Xt_„ 




3.45 




3.54 


Standard deviation of accident rate i\t,n for Poisson; ■\/At,ra(l + a\t^„) for NB), 
averaged over all values of explanatory variables ^t,n 




1.38 




2.33 


Total number of free model parameters 


22 


22 


19 


19 


Posterior average of the log-likelihood (LL) 




-2662.09:^«^«;S 




-2543.32l253S:27 


Max(LL): true maximum value of log-likelihood (LL) for MLE; maximum 
observed value of LL for Bayesian-MCMC 


-2651.16 

(true) 


-2652.37 

(observed) 


-2533.81 

(true) 


-2534.52 

(observed) 


Logarithm of marginal likelihood of data (ln[/(Y|A4)]) 




9fi79 97-2669.92 
^"'^■^'-2674.02 




-2554.16-i°:g 


Goodness-of-fit p-value 




0.000 




0.003 


Maximum of the potential scale reduction factors (PSRF) ' 




1.02304 




1.01813 


Mulli\-aria.lo potfulial sciilo rcduciiuu la(-|ur (All'SiU-') ' 




1.02 




l.()19:-!8 



Standard (conventional) Poisson and negative binomial correspondingly estimated by meiximum likelihood estimation (MLE). 

^' Standard Poisson and negative binomial correspondingly estimated by Markov Chain Monte Carlo (MCMC) simulations. 
^ The pavement quality index (PQI) is a composite measure of overall pavement quality evaluated on a to 100 scale. 

* PSRF/MPSRF are calculated separately /jointly for all continuous model parameters. PSRF and MPSRF are close to 1 for converged MCMC chains. 



Table 6.2 

Estimation results for zero-inflated and Markov switching Poisson models of annual accident frequencies 



Variable 


ZIP-x 




MSP'^ 
by MCMC 


by MLE 


by MCMC 


by MLE 


by MCMC 


/3-cocfBcients in Equation II3.4II 


Intercept (constant term) 


-6.66-?s? 

— 7.41 


-6.65-S'l7 

— 7.57 


■7 O-7.01 

— 0.54 


n oc— 6.85 

— 8.87 


-13.4-J0-7 

— lo.l 


Accident occurring on interstates 1-70 or 1-164 (dummy) 


— .595 


-.526-^f; 

— .oil 


-.594-523 

— .565 


— .586 


ao^ --544 
— .DOl -710 
— .718 


Pavement quality index (PQI) average 


-.00859-?°^ 


-.00860-o?3f 


-.0101-00822 

— .V16i) 


-.0101-00493 
— .Uioz 


nit; -.00969 


Road segment length (in miles) 


.0803 0982 


■ 0801 ?^5i 

. iU4 


.0674-0871 

.U4 1 i 


.0667-0929 

.U4Uz 


•092-ilL 

.Ud6o 


Logarithm of road segment length (in miles) 


.741-781 


■742-^89 

.Dy4 


.804-S53 


■ 808-gJ§ 

. ( 4z 


.714-776 

.6dz 


Number of ramps on the viewing side 


-.0301-00?^ 

— .U512 


-.0301-o?559 

— .U545 


-.0247-00328 

— .U4ol 


-.0247-00000181 

— .U494 


riQOO — .00825 
— .0581 


Number of ramps on the viewing side per lane per mile 


.309;^^|:5 


•308'?!; 

.2o4 


'^^"^ .233 


■ 302-2?§ 


.303 223 


Median configuration is depressed (dummy) 


■i44:gi?9 


.144-233^ 


-i49:gi5 


1 fsn-257 

■ 0443 


126-231 
' 0224 


Median barrier presence (dummy) 


-2-50:|35 


-2.5o:^:?^ 


-.0821 1 ^508 


-.828:i525 


— 9 -^0 — 1-57 


Interior shoulder presence (dummy) 










— 9 09 — 1-22 
^-^^^-2.87 


Width of the interior shoulder is less that 5 feet (dummy) 


341 ■*14 


.342:43? 


324-403 
••^^*.245 


one;. 428 
-'-'^■-'.223 


304-405 
.ciU4 204 


Outside shoulder width (in feet) 




(-vc-yo- .0358 
.U0/(5_ (jygg 


— .0635_ 0g32 


06391 :S 


- 0419--0172 

.U4iM_ gggy 


Outside barrier absence (dummy) 




17O-.0765 


_ oc:q — -1^5 
.zo.3_ 33J 


.Z0.3_ 3g5 


fjQ — . 12 1 

— 343 


Average annual daily traffic (AADT) 










-3 60^2.49 
'^-"'^-4.69 
X 10^5 


Logarithm of average annual daily traffic 


■84i:??i 


•84i:?Jt 


1 O3I-IO 
l-'^'^.959 


1.03lgll 


i-892:« 


Posted speed limit (in mph) 


ni CO. 0227 

.UIDJ QQggg 


■0164:02|io 


0082^0149 


nnsi'^70169 


•00899:017523 


Number of bridges per mile 


-.03661:0^54 




-.0246-00925 


- 0249- 00644 
.uz^a_ Q44Q 


r,noo-. 00574 
.UZZCi_ 0401 


Maximum of reciprocal values of horizontal curve radii (in 1 / mile) 


-.i48_;Jg;; 


^■148_ 201 


-.io6-;5« 


- 107^-0455 


1 r,7-.0684 
-1^'-.186 


Maximum absolute value of change in grade of a vertical curve 


.0328:011? 


.0328:0551 


.0308:0494 


0-^00-0539 


O9OS 0380 
-02U» 00341 


Number of vertical curves per roadway section 


-.14981:0^6? 


, CI -.0609 
.iOi_ 241 


,,c:7-.0362 
-11'J'-.1952 


1 1 V-.0200 
--11'-. 214 




Percentage of single unit trucks (daily average) 


.614:111 


616105 
•"1".184 


.8I43I4I5 


091 1.28 
-°^1.363 


1 001-45 
l-'-"J.548 


Number of changes per vertical profile along a roadway segment 


•0681:io|o 


.0684:ii2g 


•0398:QQf|3 


0409-0872 





Table 6.2: (Continued) 



Variable 


ZIP-x'» 


ZIP-T ^ 


MSP<^ 
by MCMC 


by MLE 


by MCMC 


by MLE 


by MCMC 


T- and -Y-coeflicients in Eauations l|3.11|l and l|3.12(l 


The model parameter r in Equation II3.11II 












Intercept (constant term) 












Logarithm of road segment length (in miles) 


- 




-1.40:1;°? 


-1 42-1-" 




Median barrier presence (dummy) 


- 






4 164-98 




Width of the interior shoulder is less that 5 feet (dummy) 


- 




- 921 -■3''° 


-1.45 


- 


Outside shoulder width (in feet) 






QQr> — . 159 

— .zzz_ 


000-. 176 




Maximum of reciprocal values of horizontal curve radii (in 1 /mile) 






.0(0 201 


.58ljfi 




Mean accident rate (At^„), averaged over all values of Xt_,i 




3.41 




3.42 


3.94 


Standard deviation of accident rate {\t,n), averaged over all 
values of explanatory variables Xt,n 




1.62 




1.67 


1.60 


Total number of free model parameters (/3-s, 7-s, a and t) 


22 


21 


25 


25 


20 


Posterior average of the log-likelihood (LL) 




-2636 m-263"-69 




— 9^1Q K/i— 2513.62 


-2149.82l2i22.28 


Max(LL): true maximum value of log-likelihood (LL) for MLE; 
maximum observed value of LL for Bayesian-MCMC 


-2625.58 

(true) 


-2626.63 

(observed) 


-2507.07 

(trLie) 


-2508.61 

(observed) 


-2080.32 

(observed) 


Logarithm of marginal likelihood of data (ln[/(Y|7V()]) 




9fi/18 /II -2643.46 




QCQQ Q-| — 252T.95 
— ZOOZ.Zi_2534 Q2 


OOOQ 07-2194.42 
ZZZM.Z( _22i4.49 


Goodness-of-fit p-valuc 




0.000 




0.009 


0.009 


Maximum of the potential scale reduction factors (PSRF) 




1.00163 




1.00252 


1.02803 


Multivariate potential scale reduction factor (MPSRF) 




1.00171 




1.00255 


1.02852 



^ Standard (conventional) ZIP-t model estimated by maximum likelihood estimation (MLE) and Markov Chain Monte Carlo (MCMC) simulations. 
^ Standard ZIP-7 model estimated by maximum likelihood estimation (MLE) and Markov Chain Monte Carlo (MCMC) simulations. 

Two-state Markov switching Poisson (MSP) model where all reported parameters are for the unsafe state s = 1. 
"1 The pavement quality index (PQI) is a composite measure of overall pavement quality evaluated on a to 100 scale. 

" PSRF/MPSRF are calculated separately/jointly for all continuous model parameters. PSRF and MPSRF are close to 1 for converged MCMC chains. 



Table 6.3 

Estimation results for zero- inflated and Markov switching negative binomial models of annual accident frequencies 



Variable 


ZINB-x'' 


ZINB-7'' 


MSNB'^ 
by MCMC 


by MLE 


by MCMC 


by MLE 


by MCMC 


f)- and a-parameters in Equation 113.711 


Intercept (constant term) 


-is.o-Jyi 

— i / .5 


-15.2-"-0 

— i / .4 


-11.6-?4^o 

— i4.» 


-ll-6"?4'6 
— 14. b 


-17.3-^?-0 


Accident occurring on interstates 1-70 or 1-164 (dummy) 


— .683 -7(1^ 


— .685 ^JJa 


7-1 c-. 602 

— .ozy 


71 C-.593 

— .000 


-.734-«i;^ 


Pavement quality index (PQI) average 


.U1ZZ_ QQggQ 


r|i or,-. 00562 


- 0140- '"'^2'' 
.U14U_ 0217 


Q221 


- 0163- °°*50 
.uloo_ 0240 


Logarithm of road segment length (in miles) 


7qi -832 
• '='^.751 


701 -829 


qoq.978 
■^^y.880 


■939:li 


007.929 
•*°'.845 


Number of ramps on the viewing side per lane per mile 


.226;?°° 


907. 306 
•^^'.149 


.298:209 


.304:394 


0-1 7.404 
••^1 ' .230 


Number of lanes on a roadway 










1 iq2-04 
1-1'^.386 


Median configuration is depressed (dummy) 




ion. 282 
■^^■^.0839 


201 -319 
•^'-'1.0820 


909.325 

.ZUZ Q7gj 




Median barrier presence (dummy) 




-1 43~1-14 






i.Da_2 45 


Width of the interior shoulder is less that 5 feet (dummy) 


QOQ.443 

.ozo 202 


noo.434 
.ozo 211 


40c. 572 
.^00 297 


407.569 
■^■^'.307 


■174.505 
••^'^.243 


Outside shoulder width (in feet) 




p.. 70-. 0207 
-■U4l'8_ 0749 


- 0532- ''-^ 
.UOOZ_ Qgg7 


— n^'^9 — -°^9 
.UOOZ_ Qgg7 


- 0537--°214 
.UOOl _.0862 


Outside barrier absence (dummy) 






-•245::^i| 


■^^'-'-.389 


- 264-- 
■^'3*-. 403 


Average annual daily traffic (AADT) 


-4.97 

X 10-5 


—4 14~3.31 
^■-^^-5.04 

X 10-5 


-1-93:6:™ 

X 10-5 


_i qi-3-16 
^•^^-5.83 

X 10-5 


-3.781^:°^ 

X 10-5 


Logarithm of average annual daily traffic 




l-yj^l.67 


i.52i:fi 


1-3^1.15 


i-95?:l| 


Number of bridges per mile 










-02141:°°^ 


Maximum of reciprocal values of horizontal curve radii (in 1/mile) 


.i4U_ 209 


_ 141 -.0734 


-.1341:0559 


-1 .JO -.0593 
217 


-io6::?if 


Percentage of single unit trucks (daily average) 


1 231-84 


-■-■^•^.646 


1 091.96 

l-'^^.693 


1 091.96 
^•'-'^.691 


1-^^.688 


Number of changes per vertical profile along a roadway segment 


nc:cc:.0930 


fic:f:2.0903 
.uouz 0226 








Over-dispersion parameter a in NB models 


•144:183 


.i50:«l 


10(1.168 
•1'^^.0925 


149. 185 
■1^^.105 


114.147 
■11^.0847 



Table 6.3: (Continued) 



Variable 


ZINB-T-'* 


ZINB-^t' 


MSNB 




by MLE 


by MCMC 


by MLE 


by MCMC 


by MCMC-: 


T- and -Y-Darameters in Eauations JS.llll and (13.1211 


The model oarameter r in Eauation l|3.11|l 


-'-•'^-2.00 


-1 70 — 1.50 
^- ''^-1.98 




- 




Intercept (constant term) 






23 l^l.S 
^'-'■^4.99 






Logarithm of road segment length (in miles) 


- 


- 


_1 34-. 942 






Median barrier presence (dummy) 


- 


- 


074.86 


4.i6i:i? 


- 


Average annual daily traffic (AADT) 






X 10"5 


10.5^^7! 

X 10^5 




Logarithm of average annual daily traffic 






o 00- 901 


-3.281-5? 




Mean accident rate {Xt,n for NB), averaged over all values of Xt,n 




3.38 




3.42 


3.88 


Standard deviation of accident rate {%/\t,n{^ + ceXt,n) for NB), 












averaged over all values of explanatory variables Xt^„ 




2.14 




2.15 


2.13 


Total number of free model parameters (/3-s, 7-s, a and t) 


16 


16 


19 


19 


16 


Posterior average of the log-likelihood (LL) 




-2510.68:2506.13 




-2436.34:2«1;12 


on/l on-2096.30 
-^i^4.»/„2153.91 


Max(LL): true maximum value of log-likelihood (LL) for MLE; 
maximum observed value of LL for Bayesian-MCMC 


-2502.67 

(true) 


-2503.21 

(observed) 


-2426.54 

(true) 


-2427.41 

(observed) 


-2049.45 

(observed) 


Logarithm of marginal likelihood of data (ln[/(Y A4)]) 




-2519.901^516.95 




-2447 S3-2443.93 


9104 91 -2186.70 
-i!l»4.^1_24g9 gg 


Goodness-of-fit p-valuc 




0.005 




0.177 


0.191 


Maximum of the potential scale reduction factors (PSRF) ° 




1.01006 




1.02200 


1.02117 


Multivariate potential scale reduction factor (MPSRF) ° 




1.01023 




1.02302 


1.02189 



^ Standard (conventional) ZINB-t model estimated by maximum likelihood estimation (MLE) and Markov Chain Monte Carlo (MCMC) simulations. 
^ Standard ZINB-7 model estimated by maximum likelihood estimation (MLE) and Markov Chain Monte Carlo (MCMC) simulations. 

Two-state Markov switching negative binomial (MSNB) model where all reported parameters are for the unsafe state s = 1. 

The pavement quality index (PQI) is a composite measure of overall pavement quality evaluated on a to 100 scale. 
" PSRF/MPSRF are calculated separately/jointly for all continuous model parameters. PSRF and MPSRF are close to 1 for converged MCMC chains. 
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Figure 6.1. The histogram of 10^ generated values for the MSNB model 
of annual accident frequencies. The vertical line shows the observed value 
of x^- 



To evaluate the goodness-of-fit for a model, we use the posterior (or MLE) es- 
timates of all continuous model parameters (/3-s, a, Po^i, Pj^q) and generate 10^ 
artificial data sets under the hypothesis that the model is truejj We find the distribu- 
tion of x^i given by equation (14.41) . and calculate the goodness-of-fit p- value for the 
observed value of x^. As a demonstration, refer to Figure 16. where the histogram 
of the generated x^ values is plotted for the MSNB model. The observed value of 
X^ is shown by the vertical line in this figure, the goodness-of-fit p-value is equal to 
the ratio of the histogram area located to the right of the vertical line and the total 
histogram area. The resulting p-values for all negative binomial models are given in 
Tables 16.11 and 16. 3[ For the ZINB-7 and MSNB models the p-values are sufficiently 
large, around 20%, which indicates that these models fit the data reasonably well. 
At the same time, for the ZINB-r model the goodness-of-fit p-value is only around 
0.5% and for the standard NB model the p-value is only around 0.3%, which indicate 
much poorer fit. Note that all Poisson models (including the MSP model) provide 



^Note that the state values S are generated by usmg pj^Xi ^^id Pi"io' where n = 1,2, N. 
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Table 6.4 

Summary statistics of explanatory variables that enter the models of an- 
nual and weekly accident frequencies 



Variable 


Mean 


Std'^ 


Min'" 


Median 


Max"* 


Accident occurring on interstates 1-70 or 1-164 (dummy) 


.155 


.363 








1.00 


Pavement quality index (PQI) average*^ 


88.6 


5.96 


69.0 


90.3 


98.5 


Road segment length (in miles) 


.886 


1.48 


.00900 


.356 


11.5 


Logarithm of road segment length (in miles) 


-.901 


1.22 


-4.71 


-1.03 


2.44 


Total number of ramps on the road viewing and opposite sides 


.725 


1.79 








16 


Number of ramps on the viewing side per lane per mile 


.138 


.408 








3.27 


Median configuration is depressed (dummy) 


.630 


.484 





1.00 


1.00 


Median barrier presence (dummy) 


.161 


.368 








1 


Interior shoulder presence (dummy) 


.928 


.258 





1 


1 


Width of the interior shoulder is less that 5 feet (dummy) 


.696 


.461 





1.00 


1.00 


Interior rumble strips presence (dummy) 


.722 


.448 





1.00 


1.00 


Width of the outside shoulder is less that 12 feet (dummy) 


.752 


.432 





1.00 


1.00 


Outside barrier absence (dummy) 


.830 


.376 





1.00 


1.00 


Average annual daily Iralhc (AADi) 


3.03 

X 10^ 


2.89 
X 10^ 


.914 

X 10^ 


1.65 
X 10^ 


14.3 

X 10^ 


Logcirithm of average annual daily traiBc 


10.0 


.623 


9.15 


9.71 


11.9 


Posted speed limit (in mph) 


63.1 


3.89 


50.0 


65.0 


65.0 


Number of bridges per mile 


1.76 


8.14 








124 


Maximum of reciprocal values of horizontal curve radii (in 1/mile) 


.650 


.632 





.589 


2.26 


Maximum of reciprocal values of vertical curve radii (in 1/mile) 


2.38 


3.59 








14.9 


Number of vertical curves per mile 


1.50 


4.03 








50.0 


Percentage of single unit trucks (daily average) 


.0859 


.0678 


.00975 


.0683 


.322 


Winter season (dummy) 


.242 


.428 








1.00 


Spring season (dummy) 


.254 


.435 








1.00 


Summer season (dummy) 


.254 


.435 








1.00 


Maximal external angle of the horizontal curve 


9.78 


12.0 





5.32 


66.7 


Outside shoulder width (in feet) 


11.3 


1.74 


6.20 


11.2 


21.8 


Number of changes per vertical profile along a roadway segment 


.522 


.908 








6.00 


Number of lanes on a roadway 


2.09 


.286 


2.00 


2.00 


3.00 


Number of ramps on the viewing side 


.310 


.865 








8.00 


Alaxiniuni a.l)S()lu[e \'alue of c-lia.nge in grade of a \'ei1ieal can\x! 


.()<)7 


1.21 








7.11 


Number of vertical curves per roadway section 


.445 


.611 








3.00 



" Standard deviation, minimum and maximum of a variable. 



^ The pavement quality index (PQI) is a measure of overall pavement quality evaluated on a to 100 scale. 
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Figure 6.2. Five-year time series of the posterior probabilities P{st^n = 
1|Y) of the unsafe state Sj_„ = 1 for four selected roadway segments 
{t = 1,2,3,4,5). These plots are for the MSNB model of annual acci- 
dent frequencies. 



relatively poor goodness-of-fit with p- value below 1% (refer to Tables WTl and [6^ . 
which can be explained by over- dispersion present in the annual frequency accident 
data. 

The estimation results also show that the over- dispersion parameter a is higher for 
the ZINB-r and ZINB-7 models, as compared to the MSNB model (refer Table [6731) . 
This suggests that over-dispersed volatility of accident frequencies, which is often 
observed in empirical data, could be in part due to the latent switching between the 
states of roadway safety. 

Now, refer to Figure 16.21 created for the case of the MSNB model (note that the 
corresponding figure for the MSP model is similar and is not reported). The four plots 
in this figure show five-year time series of the posterior probabilities P{st,n = 1|Y) 
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of the unsafe state for four selected roadway segments. These plots represent the 
following four categories of roadway segments: 

• For roadway segments from the first category we have P{st^n = 1|Y) = 1 for all 
t = 1,2,3,4,5. Thus, we can say with absolute certainty that these segments 
were always in the unsafe state St,n = 1 during the considered five-year time 
interval. A roadway segment belongs to this category if and only if it had 
at least one accident during each year (t = 1,2,3,4,5). An example of such 
roadway segment is given in the top-left plot in Figure 16.21 For this segment 
the posterior expectation of the long-term unconditional probability pi of being 
in the unsafe state is relatively large, E{pi\Y) = 0.750. 

• For roadway segments from the second category P{st.n = 1|Y) ^ 1 for all 
t = 1,2,3,4,5. Thus, we can say with high degree of certainty that these 
segments were always in the zero-accident state St^n = during the considered 
five-year time interval. A roadway segment n belongs to this category if it had 
no any accidents observed over the five-year interval despite the accident rates 
given by equation fl3.7l) were large, Xt^n ^ 1 for all t = 1, 2, 3, 4, 5. Clearly this 
segment would be unlikely to have zero accidents observed, if it were not in 
the zero-accident state all the timeo An example of such roadway segment is 
given in the top-right plot in Figure ESI For this segment E{pi\Y) = 0.260 is 
relatively small. 

• For roadway segments from the third category P{st^n = 1|Y) is neither one 



nor close to zero for all t = 1,2,3,4,5 







For these segments we cannot de- 
termine with high certainty what states these segments were in during years 



10 



Note that the zero- accident state may exist due to under-reporting of minor, low-severity accidents 



Shankar et al.l . ll997 |. 



'^-'^If there were no Markov switching, which introduces time-dependence of states via equations ()3.15|) . 
then, assuming non-informative priors 7r(st^„ = 0) = ■n{st.n = 1) = 1/2 for states st,m the posterior 
probabihties P(st.„ = 1|Y) would be either exactly equal to 1 (when At^n > 0) or necessarily below 
1/2 (when At.n = 0). In other words, we would have P{st^n — 1|Y) ^ [0.5, 1) for any t and n. Even 
with Markov switching existent, in this study we have never found any P{st^n — 1|Y) close but not 
equal to 1, refer to the top plot in Figure [ 
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t = 1,2,3,4,5. A roadway segment n belongs to this category if it had no 
any accidents observed over the considered five-year time interval and the ac- 
cident rates were not large, Xt^n ^ 1 for all t = 1,2,3,4,5. In fact, when 
Xt^n ^ 1, the posterior probabilities of the two states are close to one-half, 
P{st,n = ~ P{st^n = 0|Y) ^ 0.5, and no inference about the value of the 
state variable St^n can be made. In this case of small accident rates, the ob- 
servation of zero accidents is perfectly consistent with both states St^n = and 
St^n = 1- An example of a roadway segment from the third category is given in 
the bottom-left plot in Figure ESI For this segment E{pi\Y) = 0.496 is about 
one-half. 

• Finally, the fourth category is a mixture of the three categories described 
above. Roadway segments from this fourth category have posterior probabilities 
P{st,n = that change in time between the three possibilities given above. 
In particular, for some roadway segments we can say with high certainty that 
they changed their states in time from the zero-accident state St^n = to the 
unsafe state St^„ = 1 or vice versa. An example of a roadway segment from the 
fourth category is given in the bottom-right plot in Figure [6^ For this segment 
E{pi\Y) = 0.510 is about one-half. Thus we find a direct empirical evidence 
that some roadway segments do change their states over time. 

Next, it is useful to consider roadway segment statistics by state of roadway safety. 
Refer to Figure [631 made for the case of the MSNB model (note that the correspond- 
ing figure for the MSP model is similar and is not reported). The top plot in this 
figure shows the histogram of the posterior probabilities P{st^n = 1| Y) for all = 335 
roadway segments during all T = 5 years (1675 values of St,n in total). For example, 
we find that during five years roadway segments had P{st,n = 1|Y) = 1 and were 
unsafe in 851 cases, and they had P{st^n = 1|Y) < 0.2 and were likely to be safe in 
212 cases. The bottom plot in Figure 1^751 shows the histogram of the posterior expec- 
tations -E[p^"^|Y], where p^"^ = Pq\/ {PqXi + Pi"!o) stationary unconditional 
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Figure 6.3. Histograms of the posterior probabilities P(sj_„ = 1|Y) (the 
top plot) and of the posterior expectations i?[p^"^|Y] (the bottom plot). 
Here t = 1,2, 3, 4, 5 and n = 1,2, ... , 335. These histograms are for the 
MSNB model of annual accident frequencies. 



probabilities of the unsafe state (see Section [373]) . We find that 0.2 < |Y] < 0.8 

for all segments n = 1,2, ... , 335. This means that in the long run, all roadway 
segments have significant probabilities of visiting both the safe and the unsafe states. 



6.2 Model estimation results for weekly frequency data 

The 



12 



In this section we use weekly time periods, t = 1,2,3, ... ,T = 260 in total 
state St is the same for all roadway segments and can change every week. Four types 
of weekly accident frequency models are estimated: 



-"^^A week is from Sunday to Saturday, there are 260 full weeks in the 1995-1999 time interval. We 
also considered daily time periods and obtained qualitatively similar results (not reported here). 
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• First, we estimate the standard (single-state) Poisson and negative binomial 
(NB) models, specified by equations (13. 3p and (13. 6p . We estimate these mod- 
els, first, by the maximum likelihood estimation (MLE) and, second, by the 
Bayesian inference approach and MCMC simulations (see footnote[2]on page [60]) . 
We refer to these models as "P-by-MLE" (for the Poisson model estimated 
by MLE), "NB-by-MLE" (for NB by MLE), "P-by-MCMC" (for Poisson by 
MCMC) and "NB-by-MCMC" (for NB by MCMC). As one expects, for our 
choice of a non-informative prior distribution, the estimated P-by-MCMC and 
NB-by-MCMC models turned out to be very similar to the P-by-MLE and 
NB-by-MLE models respectively. 

• Second, we estimate a restricted two-state Markov switching Poisson model and 
a restricted two-state Markov switching negative binomial (MSNB) model. In 
these restricted switching models only the intercept in the model parameters 
vector f3 and the over-dispersion parameter a are allowed to switch between the 
two states of roadway safety. In other words, in equations (13.201) and (13.211) only 
the first components of vectors /3(q-) and /3(^) may differ, while the remaining 
components are restricted to be the same. In this case, the two states can have 
different average accident rates, given by equation (13.41) . but the rates have the 
same dependence on the explanatory variables. We refer to these models as 
"restricted MSP" and "restricted MSNB" ; they are estimated by the Bayesian- 
MCMC methods. 

• Third, we estimate a full two-state Markov switching Poisson (MSP) model and 
a full two-state Markov switching negative binomial (MSNB) model, specified 
by equations (I3.20p and (I3.2ip . In these models all estimable model parameters 
{(3-s and a) are allowed to switch between the two states of roadway safety. To 
choose the explanatory variables for the final restricted and full MSP and MSNB 
models reported here, we start with using the variables that enter the standard 
Poisson and NB models (see footnote [3] on page [601) . Then we consecutively 
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construct and use 60%, 85% and 95% Bayesian credible intervals for evaluation 
of the statistical significance of each /3-parameter. As a result, in the final 
models some components of /3(rn and are restricted to zero or restricted to 
be the same in the two stateso We do not impose any restrictions on over- 
dispersion parameters (a-s). We refer to the final full MSP and MSNB models 
as "full MSP" and "full MSNB" ; they are estimated by the Bayesian-MCMC 
methods. 

Note that the two states, and thus the MSP and MSNB models, do not have to 
exist. For example, they will not exist if all estimated model parameters turn out 
to be statistically the same in the two states, /3(o) = (which suggests the two 
states are identical and the MSP and MSNB models reduce to the standard non- 
switching Poisson and NB models respectively). Also, the two states will not exist if 
all estimated state variables St turn out to be close to zero, resulting in po-»i ^ Pi~>o 
[compare to equation (I3.23p ]. then the less frequent state St = 1 is not realized and 
the process always stays in state St = 0. 

The estimation results for all Poisson and NB models of weekly accident frequen- 
cies are given in Tables 16.51 and 16.61 respectively. Posterior (or MLE) estimates of all 
continuous model parameters (/3-s, a, Po^i and pi^o) are given together with their 
95% confidence intervals for MLE models and 95% credible intervals for Bayesian- 
MCMC models (refer to the superscript and subscript numbers adjacent to parameter 
posterior /MLE estimates in Tables ES] and ESI and see footnote O on page EI]). Ta- 
ble 16.41 on page [TO] gives summary statistics of all roadway segment characteristic 
variables Xt^„ (except the intercept). 

To visually see how the model tracks the data, consider Figure 16. 4[ The top 
plot in Figure [23] shows the weekly time series of the number of accidents on selected 
Indiana interstate segments during the 1995-1999 time interval (the horizontal dashed 
line shows the average value). This plot shows that the number of accidents per week 



^•^Of course, in the restricted models only the intercept is not restricted to be the same in the two 
states. For restrictions on other model coefBcients, see footnote H] on page [501 



77 



fluctuates strongly over time. Thus, under different conditions, roads can become 
considerably more or less safe. As a result, it is reasonable to assume that there exist 
two or more states of roadway safety. These states can help account for the existence 
of numerous unidentified and/or unobserved factors that influence roadway safety 
(unobserved heterogeneity). The bottom plot in Figure 16.41 shows corresponding 
weekly posterior probabilities P{st = 1|Y) of the less frequent state St = 1 for the 
full MSNB model. These probabilities are equal to the posterior expectations of St, 
P{st = 1|Y) = 1 X P{st = 1|Y) + X P{st = 0|Y) = E{st\Y). Weekly values of 
P{st = 1|Y) for the restricted MSNB model and for the MSP models are very similar 
to those given on the bottom plot in Figure [6^ and, as a result, are not shown on 
separate plots. Indeed, for example, the time-correlatioJ^ between P{st = 1|Y) for 
the two MSNB models (restricted and full) is about 99.5%. 

Let us now turn to model estimation results. Because estimation results for Pois- 
son models are very similar to estimation results for negative binomial models, let us 
focus on and discuss only the estimation results for negative binomial models. Our 
major findings, discussed below for negative binomial models, hold for Poisson models 
as well (unless otherwise stated). The findings are as follows. 



-"^^Here and below we calculate weighted correlation coefficients. For variable P{st — 1|Y) = i?(st |Y) 
we use weights wt inversely proportional to the posterior standard deviations of st- That is Wt oc 
min{l/std(st|Y),niedian[f/std(sf|Y)]}. 



Table 6.5 

Estimation results for Poisson models of weekly accident frequencies 



Variable 


P-by-MLE'^ 


P-by-MCMC ^ 


Restricted MSP " 


Full MSPd 


state s = 


state s = 1 


state s = 


state 3 = 1 


Intercept (c'uiisla.iil Icriiij 


-21.i_i;j 






1-21, (i 




-2o.i_-:l 


Accident occurring on interstates 1-70 or 1-164 (dummy) 


-.627I:«?« 


-.629_:7i7 


-.628::?« 






"•587I:i°? 


Pavement quality index (PQI) average*^ 


-.0132_;™f 


— 01 Q4. -^'^^^ 


- 0193- -"it^ 


~ 0193" Sil5 


- 0206 " ^52 

.u^uu_ 0252 




Road segment length (in miles) 


.0678:«»t» 


0790.0980 


0721-SS 

.u 1 ^ J- 0462 


0721 -IS J? 

.UI ^± 0462 


.UI j-i 0511 


0754 -ll??? 


Logarithm of road segment length (in miles) 


070.934 

•°' ^.810 


.862:|qq 


.862 gQ^ 


.862 go"^ 


.ouu go7 


865-oi^9 
.ouu go7 


Total number of ramps on the road viewing and opposite sides 


-.0203I:°°^f 


-. 02461 :°^^3 


- 0246~ 2i?? 


— 0246" Si?? 

■"-^^"—.0369 


- 0150-221°^ 

.uj.uu_ Q2gg 


- 0345~ 2iSS 

.uo'±u_ Q509 


Number of ramps on the viewing side per lane per mile 


oQr.471 


•402:m 


402iII 

327 


4O2HI 


•^^^^.340 


415-^?? 

■^-^"^.340 


Median configuration is depressed (dummy) 


T or". 288 
" ' .0864 


.0923 


--^■-''^.0927 


--^■^'^.0927 




349?2? 
•■^^''.180 


Median barrier presence (dummy) 


-3.051^:^? 




-3 OQ-lil 

•^-""-3.67 


-3 OO^o ii 

"^•""—3.67 


—3 11"?S? 
■-'•^^—3.78 


—3 11~?S? 
^ -3.78 


Interior shoulder presence (dummy) 


-1 11 "-445 
— 1.77 


-.980-^'^%^ 


— 2.32 


— 982-^?''.,„ 
— 2.32 


-1 12-47^*00 
-^•-^^ — 1.82 


-1 12-4T6o, 
^-^^-1.82 


Width of the interior shoulder is less that 5 feet (dummy) 


•371:271 


007. 487 
.288 


-1)07.487 
- .289 


■507.487 

..111 1 .ItMl 

.289 


.374-513 

.277 


.374-473 

.277 


Interior rumble strips presence (dummy) 


1 a'7-0''34 
— ■'-°l-.300 


_ 172-970 
' -^ — 1.30 


_ 172-967 
•l'z_j 32 


_ 172-967 
•-■-'^-1.32 






Width of the outside shoulder is less that 12 feet (dummy) 


000. 376 
•^°^.189 


979.366 
■^'^.179 


97':>.367 
-^'■^.180 


970.367 
•^'•-'.180 


27c. 369 
■^'0.185 


•276:?i 


Outside barrier absence (dummy) 


- 246 


- 2^4~'146 


.ZO'i_ 3gQ 


-^''^-.360 


-.280_;3g4 


-.280Z:llt 


Average annual daily traffic (AADT) 


-3.99:1" 
X lO""' 


X lO""' 


qr-— 3.13 
o.yo_4 g2 

X 10""' 


Q qc; — 3.13 
o.»o_4 g2 

X lO"'^ 


-3.64l2:g 

X 10""' 


-3.64Zt% 
X lO""' 


Logarithm of average annual daily traffic 




9 r)q2.27 
^•'J'^1.80 


9 (-192.26 
^•"^1.80 


9 092 -26 
^■^^1.80 


1 042. 16 
1-'^*1.73 


1 Q42.I6 

1-^*1.73 


Posted speed limit (in mph) 


01 '^l 0234 


01 zLQ-0232 
00662 


01 4Q-0232 


01 4Q-0232 
.Ui4y 00658 


09^9-9315 
.UZOZ oigg 




Number of bridges per mile 


riQiQ-. 00413 
.UZ1Z_ Q3g2 


- 0242--''°'^*'^ 


- 0243"-'^"™2 


- 0243" -""^2 


- 0254"-''°^°^ 
.UZ04_ 0427 


- 0254"-""^"7 
.UZ04_ 0427 


Maximal external angle of the horizontal curve 


nfm^fi'^-00669 


nn'^Q'S-00696 


nnQQc. 00696 
•"^"•^""-'.000917 


nO'iQ'^-00696 


00602-90822 
.UUDUZ go277 




Maximum of reciprocal values of horizontal curve radii (in 1 /mile) 


-247I:ii 


-.249Ll^f27 


-•249l117 


-.249l"|27 


— 974~-298 
■^'^-.341 


— 974 ~ -298 
•Z/4_ 34J 


Maximum of reciprocal values of vertical curve radii (in 1/mile) 


.0196:gf?i 


n-i7f;.()259 


ni7f;.0259 

.Ui (D 0Qg3Q 


ni7f;.0259 
.Ui (D 00930 


ni C9.0265 
-'Jl"^. 00998 


01 09.9265 
.UiOZ 00998 


Number of vertical curves per mile 


— .0588_;g929 


.UDZZ_ Qggg 


— Ofi9'^"-^*292 

.UDZO_ „ggg 


"•0623I:S 


-. 06441 


" -06441 


Percentage of single unit trucks (daily average) 


-'■■''^.814 


1 141-60 
J^-1^.684 


1 I4I-6O 

J^-J^^.681 


1 141-60 

1-1^.681 




1.832;g 



Table 6.5: (Continued) 



Variable 


P-by-MLE'' 


P-by-MCMC 


Restricted MSP-^ 


Full MSP'* 


state s = 


state s = 1 


state s = 


state s = 1 


Winter season (dummy) 


1 orr 254 




-.0627l1«iV3 


-.0627_i»iV3 




-364l«2'32 


Spring season (dummy) 






1 Qi .0689 
230 


101.0689 
230 






Summer season (dummy) 




-.168:!'«243 








-•345Lif68 


Mean accident rate {Xt,n), averaged over all values of Xt,„ 


— 


.0661 


.0570 


.1540 


.0533 


.1100 


Standard deviation of accident rate (At,n), averaged over all 
values of explanatory variables Xt^„ 




.1900 


.1770 


.2900 


.1730 


.2390 


Markov transition probability of jump — > 1 (po^l) 










Marktn' liaiisilioii piuijaljiiily of jiiiiii) 1 — (/>i --,o) 










Unconditional probabilities of states and 1 (po and pi) 






.902:|« and mSlUlg 


.794:871 and .206:?^ 


Total number of free model parameters (/3-s and as) 


26 


26 


27 


25 


Posterior average of the log-likelihood (LL) 




-16381.08:l«^ei:g« 


— IfiO'H^ Q7-16023.36 
^O^"^^-^' -16047.89 


-15964.02:}«^^^:^^ 


Max(LL); true maximum value of log-likelihood (LL) for MLE; 
maximum observed value of LL for Bayesian-MCMC 


-16355.68 (true) 


-16362.30 (observ.) 


-15990.70 (observed) 


-15928.03 (observed) 


Logarithm of marginal likelihood of data (ln[/(Y|jV()]) 




1000*.»(_jgggg 24 


-16056.91i1«°^°:?« 


-16001.15:;5992;86 


Goodness-of-fit p-value 




0.296 


0.404 


0.,393 


Maximum of the potential scale reduction factors (PSRF) ' 




1.02205 


1.00711 


1.00759 


Multivariate potential scale reduction factor (MPSRF) ^ 




1.02361 


1.00776 


1.00792 



* Standard (conventional) Poisson estimated by maximum likelihood estimation (MLE). 
Standard Poisson estimated by Markov Chain Monte Carlo (MCMC) simulations. 

Restricted two-state Markov switching Poisson (MSP) model with only the intercept and over-dispersion parameters allowed to vary between states. 
Full two-state Markov switching Poisson (MSP) model with all parameters allowed to vary between states. 
The pavement quality index (PQI) is a composite measure of overall pavement quality evaluated on a to 100 scale. 
' PSRF/MPSRF are calculated separately/jointly for all continuous model parameters. PSRF and MPSRF are close to 1 for converged MCMC chains. 



Table 6.6 

Estimation results for negative binomial models of weekly accident frequencies 



Variable 


NB-by-MLE* 


NB-by-MCMC ^ 


Restricted MSNB<= 


Pull MSNB <J 


state s = 


state s = 1 


state s = 


state s = 1 


Intercept (constant term) 


—91 Q — 18-7 
^^••^-23.9 




-20.9"i2'n 

— 23.0 


-19.9"oo'? 
— 22. 1 


29-7_22,8 


^^''''—22.8 


Accident occurring on interstates 1-70 or 1-164 (dummy) 


.ooo_ 






-.656_ 74g 


-.6601:568 


-.6601:568 


Pavement quality index (PQI) average 


r,iQo- 00681 
.U10Z_ 0206 


-0189I:°»t 


.U190_ Q248 


.UiyO_ 0248 


-.0220I:°l?« 


moi;--00700 

.UJ.ZO_ gig0 


Road segment length (in miles) 


■05i2:0i;i 




n^'^s-^^1-^ 

.uooo 0254 


.uooo g2g4 


■0395:0fii 


■0395:0fii 


Logarithm of road segment length (in miles) 




qnq.964 


qnn-961 
.auu g4(, 


.yuu g4o 




-9i3:ii 


Total number of ramps on the road viewing and opposite sides 


■'^-'-'^-.0327 


^Q-j — .00624 


ni 07-. 00423 
-.UlS/_ 0331 


0-107-. 00423 
-.Ui8^_ 0331 


- 


- 0264 "■""^56 
.Uz:D4_ 04S4 


Number of ramps on the viewing side per lane per mile 




■4oo:|Ii 


.397:111 


.397:111 


oc:q.429 
.oow 289 


oc:q.429 
.oo» 289 


Median configuration is depressed (dummy) 


910-314 


91 4-318 


•2ii:fJi 


•2ii:fJi 


.209:313 


.209:313 


Median barrier presence (dummy) 


O.U^_3 


_9 QQ-2-40 
z.yy_3 gy 


ni-2-42 

•^■"^-3.69 


•^•"^-3.69 


■^•"^-3.69 


•^•"^-3.69 


Interior shoulder presence (dummy) 


-'-•-'-'^-1.81 


i.UD_2 26 


-1-02^1^23 


-1-02^^23 


-i-ieii's? 


-i-ieii's? 


Width of the interior shoulder is less that 5 feet (dummy) 


Q7Q.477 
J 270 


■■J°*.279 


.386:if 


.386:if 


•380:«i 


•380:t|i 


Interior rumble strips presence (dummy) 


.iDD_ 293 


■-'-^^-1.16 


--163l«?«14 


-.163l«?«14 






Width of the outside shoulder is less that 12 feet (dummy) 


9S1 -380 


070.370 
■^'^.174 


■268:?«« 


.268:f«« 


of?7.36o 


of?7.365 
.170 


Outside barrier absence (dummy) 




.^Od_ 3gg 


.^Od_ 3gg 


.ZdO_ 3gg 


.Zdl_ 3g2 


.Zdl_ 3g2 


Average annual daily traffic (AADT) 


X IQ-s 


-4 09-^-2'' 
X 10^5 


— 4 07—3-22 

'i.U I _4 94 

X IQ-^ 


—4 07~^-22 
— 4.94 

X IQ-s 


— s qo~^-ii 

o.»o_4 72 
X 10~^ 


_4 r;o~3.61 

X 10"^ 


Logarithm of average annual daily traiJic 


2-08?:ig 


2.06?:i 


2.07?:i 


2.07?:i 




2 072-30 
^■•J' 1.84 


Posted speed limit (in mph) 


•0154:0§|t3 


ni t;n.0241 


01 fil -0251 
-'^ -'-'J -"-.00697 


01 Rl -0251 
.UlDl 00697 


ni CI .0252 
.uioi 00712 


nifil -0252 
.UlDl 00712 


Number of bridges per mile 


1-19 -, 00187 


■'^^^^-.0419 


p,r,no-. 00648 
.UZJO_ (,410 


0000-. 00648 
.UZOO_ 0410 




-.06071 :igf 


Maximum of reciprocal values of horizontal curve radii (in 1 / mile) 


-.1821:122 


■l'y-.241 


-.178_:2J9 




— .175_ 237 


— .175_ 237 


Maximum of reciprocal values of vertical curve radii (in 1/milc) 


■^''-■''-.tm;)T2 


01 77-027 
' ' ,00843 


-0183:00917 


-'^l^-^illooiV 


-0184:S, 


-oi8i:|>^Ji, 


Nuiiil>cr of \x;rti(::a.i cui\'C's pci' mile 


ossy 


.()2:i:i 
■ ''" -.0024 


.024<) 

-.().:)h(>_ 0940 


0940 




-.().-)(>.-)_:orji7 


Percentage of single unit trucks (daily average) 


1 oqI.88 
-'■■'^°.886 


-'■■25^758 


1 IqI-68 


1 IqI-68 
i.iy y0l 


.7261^28 


9 C73.39 
Z.Ofi 77 



00 



Table 6.6: (Continued) 



Variable 


NB-by-MLE 


NB-by-MCMC'^ 


Restricted MSNB "= 


Full MSNBd 


state s = 


state s = 1 


state s = 


state s = 1 


Winter season (dummy) 


i^Q.226 


■i48:iig 


-.116fs« 




.j.o»_ 269 




Spring season (dummy) 


-, 7O-.0878 
■^''^-.258 


17O-.0899 
''^-.257 


-0^32°%% 








Summer season (dummy) 


■-^ ' ^-.266 


1 on- .0963 
-.1»U_ 263 


-■0332-^146 


-.0332l"146 






Over-dispersion parameter a in NB models 


.845 


.968',% 


•S37;3g2 


1-^^.986 


■^^■J.300 




Mean accident rate (At,n for NB), averaged over all values of Xt,n 




.0663 


.0558 


.1440 


.0533 


.1130 


Standard deviation of accident rate (^At,n(l + ocXt^n) for NB), 
averaged over all values of explanatory variables Xt,n 


_ 


.2050 


.1810 


.3350 


.1760 


.2820 


Markov transition probability of jump — > 1 (po^i) 


_ 


_ 


■0933:Jgi 


•1'^°.100 


Markov transition probability of jump 1 — > (pi^o) 






•65i:!ie 


•"^'.474 


Unconditional probabilities of states and 1 (po and pi) 






.873;?29 and .127;20f3 


.798:f6i and .202;282 


Total iiuiiiljc'i' uf iitx' iiiu([t;l i)ara.iiic'l.C'is (,J-s and fi-s) 


2() 


20 


28 


28 


Posterior average of the log-likelihood (LL) 




— IfiOQ? 9-16091.3 


ic;aoi a-15807.9 


— 1 5778 (-1-1^672.9 


Max(LL): true maximum value of log-likelihood (LL) for MLE; 
maximum observed value of LL for Bayesian-MCMC 


-16081.2 (true) 


-16086.3 (observ.) 


-15786.6 (observed) 


-15744.8 (observed) 


Logarithm of marginal likelihood of data (ln[/(Y|A^)]) 




-16108.611^1°^:? 


-15850 2-i5*'*"-i 

1000U.Z_j^gg4g J 


— 1 5809 4.-15801.7 


Goodness-of-fit p-value 




0.701 


0.729 


0.647 


Maximum of the potential scale reduction factors (PSRF) * 




1.00874 


1.00754 


1.00939 


Multivariate potential scale reduction factor (MPSRF) ^ 




1.00928 


1.00925 


1.01002 



Standard (conventional) negative binomial estimated by maximum likelihood estimation (MLE). 
^ Standard negative binomial estimated by Markov Chain Monte Carlo (MCMC) simulations. 
Restricted two-state Markov switching negative binomial (MSNB) model with only the intercept and over-dispersion parameters allowed to vary between states. 



Pull two-state Markov switching negative binomial (MSNB) model with all parameters allowed to vary between states. 

® The pavement quality index (PQI) is a composite measure of overall pavement quality evaluated on a to 100 scale. 

* PSRF/MPSRF are calculated separately/jointly for all continuous model parameters. PSRF and MPSRF are close to 1 for converged MCMC chains. 
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Figure 6.4. The top plot shows the weekly accident frequencies in Indiana. 
The bottom plot shows weekly posterior probabilities P{st = 1|Y) for the 
full MSNB model of weekly accident frequencies. 



The findings show that two states exist and Markov switching models are non- 
trivial (in the sense that they do not reduce to the standard single-state models). In 
particular, we found that in the restricted MSNB model we over 99.9% confident that 
the difference in values of /3-intercept in the two states is non-zeroo In addition, 
Markov switching models (restricted and full) are strongly favored by the empirical 
data as compared to the corresponding standard models. To compare the former with 
the later, we calculate and use Bayes factors given by equation (14.31) . From Table ESI 
we see that the values of the logarithm of the marginal likelihood of the data for the 
standard NB, restricted MSNB and full MSNB models are -16108.6, -15850.2 and 



^^The difference of the intercept values is statistically non-zero despite the fact that the 95% credible 
intervals for these values overlap (see the "Intercept" line and the "Restricted MSNB" columns in 
Table l6^ . The reason is that the posterior draws of the intercepts are correlated. The statistical 
test of whether the intercept values differ, must be based on evaluation of their difference. 
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—15809.4 respectively. Thus, the restricted and full MSNB models provide consider- 
able, 258.4 and 299.2, improvements of the logarithm of the marginal likelihood as 
compared to the standard non-switching NB model. As a result, given the accident 
data, the posterior probabilities of the restricted and full MSNB models are larger 
than the probability of the standard NB model by e^^^'^ and e^^^'^ respectively!^ Note 
that we use equation (14. 2 p for calculation of the values and the 95% confidence inter- 
vals of the logarithms of the marginal likelihoods reported in Tables l675l and WM The 
confidence intervals are found by bootstrap simulations (see footnote [7] on page [621) . 
We can also use a classical statistics approach for model comparison, based on the 



maximum likelihood estimation (MLE). Referring to Table [631 the MLE gives the 
maximum log-likelihood value —16081.2 for the standard NB model. The maximum 
log-likelihood values observed during our MCMC simulations for the restricted and 
full MSNB models are —15786.6 and —15744.8 respectively. An imaginary MLE, at 
its convergence, would give MSNB log-likelihood values that would be even larger 
than these observed values. Therefore, if estimated by the MLE, the MSNB models 
would provide very large (at least 294.6 and 336.4) improvements in the maximum 
log-likelihood value over the standard NB model. These improvements would come 
with only modest increases in the number of free continuous model parameters [jS-s 
and a-s) that enter the likelihood function. Both the Akaike Information Criterion 
(AIC) and the Bayesian Information Criterion (BIC) would strongly favor the MSNB 
models over the NB model (see footnote [8| on page [62!) . 

To evaluate the goodness-of-fit for a model, we use the posterior (or MLE) es- 
timates of all continuous model parameters (/?-s, a, Po^i, P\^o) and generate 10^ 
artificial data sets under the hypothesis that the model is trucj. We find the distri- 
bution of x^, given by equation (14. 4p . and calculate the goodness-of-fit p- value for the 



^^In addition, we find DIG (deviance information criterion) values 32219, 31662, 31577 for the NB, 
restricted MSNB and full MSNB models respectively. We also find DIG values 32771, 32086, 31946 
or the Poisson, restricted MSP and full MSP models respectively. This means that the MSNB 
(MSP) models are favored over the standard NB (Poisson) model [the full MSNB (MSP) is favored 
most]. However, we prefer to rely on the Bayes factor approach instead of the DIG (see footnote [51 
on page lSTj) . 

^^Note that the state values S are generated by using po^i and Pi^o- 
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observed value of x^- The resulting p- values for the NB models are given in Table \QM 
These p-values are around 65-70%. Therefore, all models fit the data well. 

Focusing on the full MSNB model, which is statistically superior because it has 
the maximal marginal likelihood of the data, its estimation results show that the less 
frequent state = 1 is about four times as rare as the more frequent state St = 
[refer to the estimated values of the unconditional probabilities pQ and pi of the states 
and 1, which are given by equation (13.161) and reported in the "Full MSNB" columns 
in Table 16. 6j . 

Also, the findings show that the less frequent state Sf = 1 is considerably less safe 
than the more frequent state St = 0. This result follows from the values of the mean 
weekly accident rate Xt^n [given by equation (13.71) with model parameters /3-s set to 
their posterior means in the two states], averaged over all values of the explanatory 
variables X( „ observed in the data sample (see "mean accident rate" in Table . 
For the full MSNB model, on average, state St = 1 has about two times more accidents 
per week than state St = haso Therefore, it is not a surprise, that in Figure 16.41 
the weekly number of accidents (shown on the top plot) is larger when the posterior 
probability P{st = 1|Y) of the state St = 1 (shown on the bottom plot) is higher. 

Note that the long-term unconditional expectation of accident frequency At^n is 
E{At^n) = Po(Ai°i)^ +Pi(a2)^, where aJ°2 = exp(/3(o)Xi,„) and aJ;^ = exp(/3(i)Xi,„) 
are the mean accident rates in the states St = and St = I respectively [see equa- 
tion (13. 7p ]. and {■ ■ ■)t means averaging over time. The unconditional expectation 
E{At^n) should be used in all predictions of long-term averaged accident rates on the 
n^^ roadway segment. In the formula for this expectation, the mean accident rate A^^^ 
is averaged over the two states by using the stationary unconditional probabilities po 
and pi (see the "unconditional probabilities of states and 1" in Table [631) . 



^^Note that accident frequency rates can easily be converted from one time period to another (for 
example, weekly rates can be converted to annual rates). Because accident events are independent, 
the conversion is done by a summation of moment-generating (or characteristic) functions. The sum 
of Poisson variates is Poisson. The sum of NB variates is also NB if all explanatory variables do not 
depend on time (Xt_„ — X„). 
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It is also noteworthy that the number of accidents is more volatile in the less 
frequent and less-safe state (st = 1). This is reflected in the fact that the standard 



deviation of the accident rate (stdj^„ = At,„(l + aA(,„) for NB distribution), av- 
eraged over all values of explanatory variables Xt^„, is higher in state = 1 than 
in state = (refer to Table [631) . Moreover, for the full MSNB model the over- 
dispersion parameter a is higher in state St = \ [a = 0.443 in state = and 
a = 1.16 in state st = 1). Because state = 1 is relatively rare, this suggests that 
over-dispersed volatility of accident frequencies, which is often observed in empirical 
data, could be in part due to the latent switching between the states, and in part due 
to high accident volatility in the less frequent and less safe state St = 1. 

To study the effect of weather (which is usually unobserved heterogeneity in most 
data bases) on states. Table 16.71 gives time-correlation coefficients between poste- 
rior probabilities P{st = 1|Y) for the full MSNB model and weather-condition vari- 
ables. These correlations were found by using daily and hourly historical weather 
data in Indiana, available at the Indiana State Climate Office at Purdue University 
(www.agry.purdue.edu/climate). For these correlations, the precipitation and snow- 
fall amounts are daily amounts in inches averaged over the week and across several 
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The 



weather observation stations that are located close to the roadway segments 
temperature variable is the mean daily air temperature (°F) averaged over the week 
and across the weather stations. The effect of fog/frost is captured by a dummy 
variable that is equal to one if and only if the difference between air and dewpoint 
temperatures does not exceed 5°F (in this case frost can form if the dewpoint is be- 
low the freezing point 32° F, and fog can form otherwise). The fog/frost dummies 
are calculated for every hour and are averaged over the week and across the weather 
stations. Finally, visibility distance variable is the harmonic mean of hourly visibility 



Snowfall and precipitation amounts are weakly related with each other because snow density 
(g/cm^) can vary by more than a factor of ten. 
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Table 6.7 

Correlations of the posterior probabilities P{st = 1|Y) with weather- 
condition variables for the full MSNB model 





All year 


Winter 

(Nov.-Mar.) 


Summer 

(May-Sept.) 


Precipitation (inch) 


0.031 




0.144 


Temperature {°F) 


-0.518 


-0.591 


0.201 


Snowfall (inch) 

> 0.2 (dummy) 


0.602 


0.577 




0.651 


0.638 




Fog / Frost (dummy) 


0.223 


(frost) 0.539 


(fog) 0.051 


Visibility distance (mile) 


-0.221 


-0.232 


-0.126 



distances, which are measured in miles every hour and are averaged over the week 
and across the weather stations]^ 

Table W77\ shows that the less frequent and less safe state = 1 is positively corre- 
lated with extreme temperatures (low during winter and high during summer), rain 
precipitations and snowfalls, fogs and frosts, low visibility distances. It is reasonable 
to expect that during bad weather, roads can become significantly less safe, resulting 
in a change of the state of roadway safety. As a useful test of the switching between 
the two states, all weather variables, listed in Table I6.7[ were added into our full 
MSNB model. However, when doing this, the two states did not disappear and the 
posterior probabilities P{st = 1|Y) did not changed substantially (the correlation 
between the new and the old probabilities was around 90%). As another test, we 
modified the standard single-state NB model by adding the weather variables into 
it. As a result, the marginal likelihood for this model improved noticeably, but the 
modified single-state NB model was still strongly disfavored by the data as compared 
to the restricted and full MSNB models. This result emphasizes the importance of 
the two-state approach. 

^°The harmonic mean Jof distances dn is calculated as ~ (V^) 12n=i '^n^^ assuming c?„ = 0.25 
miles if (i„ < 0.25 miles. 
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Let us give a brief summary of the effects of explanatory variables on accident 
rates. We will focus on those variables that are significantly different between the 
two states in the full MSNB model. Table 16.61 shows that parameter estimates for 
pavement quality index, total number of ramps on the road viewing and opposite 
sides, average annual daily traffic (AADT), number of bridges per mile, percentage of 
single unit trucks, and season dummy variables are all significantly different between 
the two states. All these differences are reasonable and could be explained by ad- 
verse weather/pavement conditions in the less-safe state st = I, and by the resulting 
lighter-than-usual traffic and more alert /defensive driving in this state. In particular, 
as compared to variable effects in the safe state St = 0, in the less safe state St = 1 
an improvement of pavement quality leads to a smaller reduction of the accident 
rate, an increase in percentage of single unit trucks results in a larger increase of the 
accident rate, and an increase in AADT leads to a smaller increase of the accident 
rate (note that the effects of AADT and its logarithm should be considered simul- 
taneously). An increase in number of ramps and bridges, and the summer season 
indicator significantly reduce the accident rate only in the less-safe state St = 1. The 
winter season indicator reduces the accident rate only in the safe state St = (this 
result, which might look counter- intuitive, could be explained by an increase in cases 
of over-confident, reckless driving during good weather/pavement conditions, unless 
there is a winter). 

Finally, because the time series in Figure 16.41 seem to exhibit a seasonal pattern 
[roads appear to be less safe and P{st = 1|Y) appears to be higher during winters], we 
estimated MSNB and MSP models in which the transition probabilities po-,i and pi^o 
are not constant (allowing each of them to assume two different values: one during 
winters and the other during non- winter seasons) o However, these models did not 



^^Let us briefly describe how these models can be specified by using the general representation of 
Markov switching models, given in Section [5.21 We define the winter seasons to be from November 
to March. The non-winter seasons are from April to October. For relations between the real time 
indexing and the auxiliary time indexing we have i = t, T = T, n = n, — N, T_ = {}. The 
elements of set T = {1, 14, 45, 67, 97, 119, 149, 171, 201, 223, 254, 261} are in weekly time units and 
contain the left boundaries of the winter and non-winter time intervals for the years 1995-1999. The 
total number of time intervals is R — 11. Transition probabilities Po^i, Pi'l^o, Po'^i ^^'^ Pi^c 
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perform as well as the MSNB and MSP models with constant transition probabilities 
[as judged by the Bayes factors, see equation (14.31) ] J^^ 



which are for the first winter and first non-winter intervals are free parameters. All other transition 
probabilities are not free: for the remaining winter intervals they are restricted to Pq^j^ and p^^q, 

(2) (2) 

and for the remaining non-winter intervals they are restricted to Pq_^i and Pi_,o- 

^^We have only six (five full) winter periods in our five-year data. MSNB and MSP with seasonally 

changing transition probabilities could perform better for an accident data that covers a longer time 

period. 
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CHAPTER 7. SEVERITY MODEL ESTIMATION RESULTS 



In this chapter we present model estimation results for accident severities. We esti- 
mate a standard multinomial logit (ML) model and a Markov switching multinomial 
logit (MSML) model. We compare the performance of these models in fitting the 
accident severity data. 

The severity outcome of an accident is determined by the injury level sustained 
by the most injured individual (if any) involved into the accident. In this study we 
consider three accident severity outcomes: "fatality", "injury" and "PDO (property 
damage only)", which we number as i = 1,2,3 respectively (/ = 3). We use data 
from 811720 accidents that were observed in Indiana in 2003-2006, and we use weekly 
time periods, t = 1, 2, 3, . . . , T = 208 in totalQ The state St can change every week. 
To increase the predictive power of our models, we consider accidents separately 
for each combination of accident type (1-vehicle and 2- vehicle) and roadway class 
(interstate highways, US routes, state routes, county roads, streets). We do not 
consider accidents with more than two vehicles involved!^ Thus, in total, there are 
ten roadway-class-accident-type combinations that we consider. For each roadway- 
class-accident-type combination the following two types of accident frequency models 
are estimated: 

• First, we estimate a standard single-state multinomial logit (ML) model, which 

is specified by equations fl3.13p and (13.141) . We estimate this model, first, by the 

maximum likelihood estimation (MLE), and, second, by the Bayesian inference 

approach and MCMC simulations [for details on MLE modeling of accident 

week is from Sunday to Saturday, there are 208 full weeks in the 2003-2006 time interval. 
^Among 811720 accidents 241011 (29.7%) are 1-vehicle, 525035 (64.7%) are 2-vehicle, and only 45674 
(5.6%) are accidents with more than two vehicles involved. 
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severities see |Malyshkinal . l2006l |: see also footnote [2] on page [60]. We refer to 
this model as "ML-by-MLE" if estimated by MLE, and as "ML-by-MCMC" if 
estimated by MCMC. As one expects, for our choice of a non- informative prior 
distribution, the estimated ML-by-MCMC model turned out to be very similar 
to the corresponding ML-by-MLE model (estimated for the same roadway-class- 
accident-type combination). 

• Second, we estimate a two-state Markov switching multinomial logit (MSML) 
model, specified by equation (13.241) . by the Bayesian-MCMC methods. To 
choose the explanatory variables for the MSML model, we start with using 
the variables that enter the standard ML model (see footnote [3] on page [60l) . 
Then, we consecutively construct and use 60%, 85% and 95% Bayesian credible 
intervals for evaluation of the statistical significance of each /5-parameter. As 
a result, in the final model some components of /3(o) and (3f^i-^ are restricted to 
zero or restricted to be the same in the two states (see footnote [Hon pageEHD- 
We refer to this model as "MSML" . 

Note that the two states, and thus the MSML models, do not have to exist for 
every roadway-class-accident-type combination. For example, they will not exist if 
all estimated model parameters turn out to be statistically the same in the two states, 
/3(o) = (which suggests the two states are identical and the MSML models reduce 
to the corresponding standard ML models). Also, the two states will not exist if all 
estimated state variables St turn out to be close to zero, resulting in po^i ^ Pi~^o, 
compare to equation (13.261) . then the less frequent state = 1 is not realized and the 
process stays in state St = 0. 

Turning to the estimation results, our findings show that two states of roadway 
safety and the appropriate MSML models exist for severity outcomes of 1-vehicle ac- 
cidents occurring on all roadway classes (interstate highways, US routes, state routes, 
county roads, streets), and for severity outcomes of 2- vehicle accidents occurring on 
streets. The model estimation results for these roadway-class-accident-type combina- 



91 



tions, where Markov switching across two states exists, are given in Tables I7.H - I7.6[ 
We do not find existence of two states of roadway safety in the cases of 2-vehicle 
accidents on interstate highways, US routes, state routes and county roads (in these 
cases all estimated state variables St were found to be close to zero, and, therefore, 
MSML models reduced to standard non-switching ML models). The standard ML 
models estimated for these roadway-class-accident-type combinations are given in Ta- 
bles IA.lHA.4l in the Appendix. In Tables I7.1H7.6I and Tables IA.lHA.4l posterior (or 
MLE) estimates of all continuous model parameters (/?-s, po-+i and pi^o) are given 
together with their 95% confidence intervals (if MLE) or 95% credible intervals (if 
Bayesian-MCMC), refer to the superscript and subscript numbers adjacent to param- 
eter posterior/MLE estimates, and also see footnote [5] on page |6T1 Table 17.71 gives 
description and summary statistics of all accident characteristic variables ,j except 
the intercept. 

Because we are mostly interested in MSML models, below let us focus on and 
discuss only model estimation results for roadway-class-accident-type combinations 
that exhibit existence of two states of roadway safety. These roadway-class-accident- 
type combinations (six combinations in total) include cases of 1-vehicle accidents 
occurring on interstate highways, US routes, state routes, county roads, streets, and 
2-vehicle accidents occurring on streets, see Tables I7.1H7.6I 

The top, middle and bottom plots in Figure mi show weekly posterior probabilities 
P{st = 1|Y) of the less frequent state St = 1 for the MSML models estimated for 
severity of 1-vehicle accidents occurring on interstate highways, US routes and state 
routes respectively^ The top, middle and bottom plots in Figure 17.21 show weekly 
posterior probabilities P{st = 1|Y) of the less frequent state St = 1 for the MSML 
models estimated for severity of 1-vehicle accidents occurring on county roads, streets 
and for 2-vehicle accidents occurring on streets respectively. 



■^Note that these posterior probabihties are equal to the posterior expectations of St, P{st = 1|Y) = 
1 X P{st = 1|Y) + X Pist = 0|Y) = E{st\Y). 
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Table 7.1 

Estimation results for multinomial logit models of severity outcomes of 
one-vehicle accidents on Indiana interstate highways 



Variable 


ML-by-MLE 


ML-by-MCMC 


MSML 


state 3 = 


state 3 = 1 


fatality 


injury 


fatality 


injury 


fatality 


injury 


fatality 


injury 


intercept 


-n q-i° i 


-s.mzlil 




70—3.56 
'^-3.88 


19 9-10.5 

1^-^-14.4 


-3-98l4:I? 


19 9-10.5 

-l^-^_14.4 


Q QQ — 2.98 
— O.ZZ_3 45 


sum 


o^c.329 
.zoo ]^42 


9QC.329 
.zoo 142 


007.329 
.zor 143 


007.329 
•^"^'.143 


-176f,ii 


176-293 
•1'".0551 


-i76:?,ii 


.615;232 


thday 


■ '^°-1.48 


- 


-853:i^"5« 


- 


-.8721,225 


- 


-.8721,225 


- 


cons 


■^''°-.623 


— 418~-213 




-4251:^^^ 


-.5661 


-.5661 


-.5661 




light 


QQf) — .0368 
— 74g 


1 q>7.224 

•-'-'^'.0501 


QQ-7-.0301 

-.3»C_ 740 


140.230 

•l-*'^.0568 


o'jQ — -0236 


1 QQ.226 
•l-'^^.0522 


070 — .0236 
.0(0_ 729 


1 QQ.226 
•1-'^^.0B22 


precip 


1 QQ — .830 
-i.d8_j^ 92 


-36iiir7 


1-^1-1.99 


-.3631:26^ 


-1-541^:?? 


.OUO_ 729 


i.04_2 10 




slush 


-i.^»_2 4g 


-•432l;i^ 


-1-431214 


-■438Z:lg 


-. 05151 :f7l 


-.0515lj6i 


-.05i5::g« 




driv 


cy-i .929 
.0^1 213 




-577:ii 




.566:ifi' 




.566:1?? 




curve 


•ii4:g;is 


114-212 


116-213 
■11°. 0186 


116-213 

■110. 0186 










driver 


4.24i:f0 


J.. OOj 43 


4-39i:i| 


■'■•^*1.43 


4-48i:Ii 


2.00?i| 


4-48i:Ii 


7-1 c. 946 
•'■1''.468 


hl20 


■ '^'^.693 


700.887 
- '^^.693 


-790:111 


790-891 
•'^^.691 


70c:. 886 


70 886 
-'=5.684 


7Qt:.886 


70t:.886 


moto 


3-883:17 


2 743.12 
^- '^2.36 


3.87t:g 


2-752:37 


4 61 ■"''•■la 
4.013 74 


q 9-33.83 
O.ZO2 70 




1 •?q2.49 
l-'^y.326 


vage 


.0285;S|2g5 


.0285;{|2™ 


.0286;Jj|^'j' 


.0286:«3™ 




.0286;[|']oJ 




.0286;o2oo 


X27 


■366:ti 


1 90.159 


-367;it 


•1^'^.0861 


■366:i| 


1 9/1.161 
■1^^.0874 


■366:i| 


1 9/1.161 
■1^^.0874 


rmd2 


2.60f:00 




2-86f:i 




2.86f:ii 




2.86f:|i 




^33 


1.24212 


-•^^^-.665 


1 iq2.02 
J^--'-°.206 


- Q4t;--033B 


l.DO g2i 


QQO-.0198 
.OOZ_ ggg 




009-. 0198 

.OOZ_ ggg 


^35 




000.410 

.ozo 246 




001 .413 
.001 24S 




994-338 
-^^^.107 




47q.637 
-^'^.328 








.00724 


.176 


.00733 


.174 


.00672 


.192 


PO^l 






1 CI .254 
•1°1.0704 








000.632 
.OOU ig4 


/;() and pi 






.(iS:-!;^:i,; and .:-!17;|™ 


# free par. 


25 


25 


28 


averaged LL 




-8486.78:;:g2:«? 


-8396.78lg™:2i 


max{LL) 


-8465.79 (true) 


-8476.37 (observed) 


-8358.97 (observed) 


marginal LL 




-8498.46I««^:i 


-8437.07_g4|Q;o2 


Good.-of-fit 




0.255 


0.222 


max(PSRF) 




1.00302 


1.00060 


MPSRF 




1.00325 


1.00067 


# observ. 


accidents = fatalities + injuries + PDOs; 19094 = 143 + 3369 + 15582 
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Table 7.2 

Estimation results for multinomial logit models of severity outcomes of 
one-vehicle accidents on Indiana US routes 



Variable 


ML-by-MLE 


ML-by-MCMC 


MSML 


state s = 


state s = 1 


fatality 


injury 


fatality 


injury 


fatality 


injury 


fatality 


injury 




,. .",.()() 






-2.1241:: 




--) n- -L-7i 


-'■).72 


-2.79__rj::]J 


sum 


■514-?« 


■200-S!27 


.509-?83 


.200-305 
.uyoi 


• 190 n*7on 

.u f oy 


• 190 ri'Ton 


• 190 rk^on 
.U * oy 




light 


-498lii 


1 04-287 


-■492I:^g 


■203:??« 


-.4931:130 


107-290 




107.290 

■1'^'.105 


snow 






l.ou_2 47 




-1.10:2^^1 


■-'-"'^.0115 


-1-10 = 2'27 


165-317 
■l'3°.0115 


nojun 


■ ''-'^.149 


oi '7.335 
■^^'.0994 


797I.3I 

■ .199 


91 0.331 
■^^•-'.0968 


•787i2i 


91 4-332 
■^J^^.0965 


•787121 


91 4.332 
•^^^.0965 


str 




_ oqtr-.lQl 


— 739" -377 


— .296_ ggg 


0*7 — .372 
'••^'-1.09 


- 294-189 


-7-37:i3^^ 


— 994-18^ 
■^^*-.398 


env 




-i-89il:^i 


— ^ "^1 ~2.81 
'J-^'-'-_4.32 


1.09_2 00 


_o r;Q-2.89 
•J.J»_4.40 


-2-09I^:8| 


o.o»_4 40 


-.701I^?3 


hllO 




504-68! 


.562:1^0 


.562:650 


.560:048 


.560:048 


.560:048 


.560:048 


moto 


2-62f:fI 


0.ZU2 


2-57?:ii 


-391 3.56 
o.z±2.87 


Q 993.58 
•^■^^2.88 


3-222:88 


993.58 
•-'■^^2.88 


q 993.58 


vage 


.UJDJ 0283 


r)oco.()444 


.UOOl 0287 


■'J'-'°'.0287 




.U,50D 0285 




•0366:«4|7 


X29 


.UODO Qoggo 




nQ7Q-0643 
.UOlO Q]^]^7 


ni 1 O.0176 
.UllO Q06I6 


■0285:ofo4 


ni 02-0178 
.UiU/ 00635 




.uizu 00635 


r21 






000.0517 

.Z/;0_ ggg 


O9Q.0517 

.^ZO_ ggg 


O9A.0504 


004.0604 


004 0504 
"•^^^-.401 


OO4.0504 


^33 






1 

J^-1'^.315 




-1 97I.98 
-■-■^'.452 




1 071-98 
^•^'.452 






01 1 40213 




Qi 1.^.0211 




nini.0200 
.UiUi 0000542 








wday 




^.104fi« 




_ 104-0124 




1 9c:. 0242 
.±zo_ 227 










079. 362 
■^'^.183 




■^'".186 




.280:309 




.280:309 








.00717 


.171) 


,00iS2:i 


.18:1 


.00218 


.158 








.0767:ii^g 


Pl->0 






■613:1^ 


po and pi 






.887:^70 and -113:0409 


# free par. 


24 


24 


25 


averaged LL 




-7406.39I?i;^:«i 


-7349.061^335.46 


max(LL) 


-7384.05 (true) 


-7396.37 (observed) 


-7318.21 (observed) 


marginal LL 




-7417 q8-''*l3.72 
— ^41C.yO_742() 23 


—TVTl 4Q— 7369.62 
'•J' ' -^^-7380.00 


Good.-of-fit 




0.337 


0.255 


max(PSRF) 




1.00319 


1.00073 


MPSRF 




1.00376 


1.00085 


# observ. 


accidents = fatalities + injuries + PDOs: 17797 = 138 + 3184 + 14485 
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Table 7.3 

Estimation results for multinomial logit models of severity outcomes of 
one-vehicle accidents on Indiana state routes 



Variable 


ML-by-MLE 


ML-by-MCMC 


MSML 


state s = 


state 3 = 1 


fatality 


injury 


fatality 


injury 


fatality 


injury 


fatality 


injury 


intercept 




^•°'-1.80 


, „„ Q 71 

-4-03_3:^i 


— 1 71 


„ . . Q 1 n 

'^•**-3.79 


-i.68iJ:|t 


-^■96i^:i| 


-i.68i}:|4 


sum 




.232;f|-;(^ 


.232:f07 


.232:307 


■238:3i| 


^ ^ ^ Q 1 /I 

■^•^".163 


■238:3i| 


•238:314 


Xrz 




'MY? 

-.390_:f^^ 


-.395_:30« 


-.395_:3g!; 




00c — .296 

.OO0_ 474 


— 2.U5_3 g2 


q Qc; ..^.-"u 
"J*0'J__,474 


light 


-.646_ 


.i93:?« 


-64i::«^ 


.199:^7 


-.689I:^^« 




-689I:*|8 


977.378 
'.177 


prccip 


-•854f«4 




-.8681,^1^ 




24 




-.8291,4^8 




driv 


r 00 — .225 




-.5961:^50 




- 58q--241 




•Ooy_.960 




str 


-2841:^1^ 


-284::^^^ 


-.2831:^1^ 


-■283I:ii^ 


iiy-.0184 

■-^-^'-.214 


11 7-. 0184 
-.214 


11 7-. 0184 
-.214 


-4651:3^0 


env 


-4.23-i« 


-1.83-:™ 


-4.28:3:67 


-l-84li>?^ 


-4 40-3.79 


-2.3o:^:« 


-4.401-™ 


-i.4iii:S 


hl20 


.840:?g 


•840:?g 


.863:^1 J 


.863:^1 J 




■86i:?f| 


1 -64:^811 


.86i:?44 




in3.3i 
■'•-^'^2.89 


3 in^-si 

•^■^"2.89 


inS.Sl 
•^■^"2.89 


in3.31 
•^■^"2.89 


q 073.66 
■^••^'3.09 


q 073.66 
'^•'^'3.09 


q 073.66 
'-*-'^'3.09 


2.823:19 


X.27 


•0557:S 


.uoo ( 0266 


.uooo 027g 


n^fi^-0858 


■0942:1383 


■0942:1383 


■0942:1383 










i.87?:i 


447.768 

.124 


i-87?:« 


461-782 


i-87?:« 


•461:11? 


^34 


14 fi21-4 

X 10^3 


-2.80Z4«",JJ 
X 10-3 


X 10-3 


—9 71 —-723 
^•'^-4.69 

X 10-3 


14.5|i6l 

X 10-3 


^•^0-4.44 
X 10-3 


14.5|i6| 

X 10-3 


-2 46- 'l'^3 
Z-^D_4 44 

X 10-3 


X35 


— 4Qfi^''^ll 
•^^"-.780 


070.344 
.z/y 214 




070.343 
•^'°.213 


_ 473 --192 
.'±(0_ 7g4 


nQq..348 
.zoo 218 


_ 470-.I92 
■^''^-.764 


OOQ.348 
.zoo 218 


vagc 




0004.0392 




.0335:q27 7 




nq'30.0390 
.uooz Q274 




0009.0390 
.uooz Q274 


othUS 








•^^^-.679 




-.4361:208 




-.4361:208 








.0089 


.179 


.00951 


.180 


.00804 


.179 


pa- -i 














•450:ij;| 


pQ and pi 






.574:|8i and .426:|?3 


# free par. 


22 


22 


28 


averaged LL 




-13867.40li3|6i.92 


iq7ft1 7^ — 13765.02 
J-0(Ol. /O_i3goo gg 


rnax(LL) 


-13846.60 (true) 


-13858.00 (observed) 


-13745.61 (observed) 


marginal LL 




-13877.89_J3||o:38 


-13820.20i13|°?:85 


Good.-of-fit 




0.515 


0.445 


max(PSRF) 




1.00027 


1.00029 


MPSRF 




1.00041 


1.00045 


# observ. 


accidents = fatalities + injuries + PDOs: 33528 = 302 + 6018 + 27208 
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Table 7.4 

Estimation results for multinomial logit models of severity outcomes of 
one- vehicle accidents on Indiana county roads 



Variable 


ML-by-MLE 


ML-by-MCMC 


MSML 


state s = 


state s = 1 


fatality 


injury 


fatality 


injury 


fatality 


injury 


fatality 


injury 


intercept 




-^■""^-1.71 


"•^'^-7.12 


— 1 65-^^^® 


-6.161^:?^ 


-i.8i_l;93 


— 7 51 -6.75 


1Q-1.99 
^.1^0-2.26 


sum 


1 "il -201 


1 "il -201 


14Q.2OO 
•-^*^.0988 


14Q.2OO 
■^^'^.0988 


140.194 

■^^''.0891 


140.194 

■^^^.0891 


— 


140.194 
•-^*^.0891 


wday 




-. 09871 


07c: — -102 
.^'0_.446 


-.09521;°^°'^ 


- 14fi--''934 

■-'-^'J-.igs 


— 14fi--''^34 
-.198 


-.i46i:?|f 


— 


dayt 


-•456::^g 




-■4431:^^ 




■^''^-.709 


— 




— 


X\2 






-.667:i2?s' 


-.169I:«^f 


-.689:i227 


-207::«^(5i 


-■68911^^7 


— 


slush 




-■2931 


^■^^^^-1.68 


-•294:Jg 


-97811^:^1 


-■2901:21? 


-97811^:^1 


-.2901:^12 


nojun 


.4i8:f|i 




.165 




.203 


267-331 
.203 


267 331 
.203 


267331 
■-'"'.203 


env 




-1 40-1-^^ 


Q — ■ 3 . 23 

— a. ( 1_4 25 


^^■''U.i 45 


-3.711^:23 


-1 76-1-*'^ 

^•'"-1.84 


-3.711^:1 


700 — .634 
■ '•^'^-.830 


hl20 


-'■•'^^1.08 


■825; 


^■■5^1.10 


■814.767 


1.34l;fg 


-809;?i5 


1.34l;f9 


■809;?i 


moto 


3-03i:i 


2 7q2-95 


•^■"^2.66 


^■'"2.62 


2.893;05 


2-89i;?i 




2.893:"^ 


vagc 




.IMDU (1322 


■"-'-'",()()27(i 


Qof;i .0398 
■"■'"-'-, 0323 


-"-^•'■^,()()1()4 


noro.0391 
■''■'■'■', ()31(i 


■"-'-•'■^,(}()1()4 


.0353;»3fi 


A 27 


.104 




■"'07 - i'' 


■-'--'-".U947 


.20();j;i^ 


--'--'-".l)9(,(> 


■2oi);i:i4 






nioc.0279 

.UlOd QQgiQ 




ni sc. 0280 
.Ui»b 00927 




ni s'? 0278 
.UlOO 00901 




nis'?.0278 
.UlOO 00901 




V" 

-^33 


o 9n2.57 
^•^^1.86 


.748|t? 


01 2.56 
^■^-■-1.84 


740.942 
•'*"^.543 


2.342;<Ii 


71 e. 916 
■ '10.516 




716-916 

■'1".516 


^34 


13.4^«iJ 
X 10"^ 


X 10"^ 


-■-3.48 07 

X IQ-^ 


-5 56-''-12 

O.OD_7 00 

X 10"'^ 


9.99^^96 
X 10"'^ 


— 5 1 7-3.73 
■^■1 ' -6.62 

X 10"'^ 


o 11 4.45 
3. 111. 73 

X 10"^ 


_P. ,7-3.73 
^■-■-'-6.62 

X 10"-^ 


^35 




■^^'^.203 


-■3621:";; 


.248:i;.i 


-''°^-.581 


.220;?^! 


- 384"-"^2 
■■^"^-.581 


.319;«? 


day 




inc. 147 
•^"°.0626 




.124;""' 0813 




■108;o650 




■108;o650 


str 




_ 147-.IOI 
194 




-i46i:?^r 




-.o8ioi:?3g'' 












.00945 


.227 


.0102 


.226 


.00594 


.228 








■0780;i3|g 








.324;«i 


po and pi 






.803;|02 and .197;3262 


# free par. 


30 


30 


34 


averaged LL 




— 30740 2Q-'^0733.70 
'^"'*"-''^-30748.77 


-30513 QS-30499.38 
ouoio.ao_3053o.oo 


max(LL) 


-30666.16 (true) 


-30728.43 (observed) 


-30480.05 (observed) 


marginal LL 




-30754 24-30^''9-02 


— 30547 OQ-30535.46 


Good.-of-fit 




0.242 


0.303 


max(PSRF) 




1.00080 


1.00025 


MPSRF 




1.00098 


1.00041 


# observ. 


accidents = fatalities + injuries + PDOs: 60782 = 581 + 13797 + 46404 
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Table 7.5 

Estimation results for multinomial logit models of severity outcomes of 

one- vehicle accidents on Indiana streets 



Variable 


ML-by-MLE 


ML-by-MCMC 


MSML 


state s = 


state s = 1 


fatality 


injury 


fatality 


injury 


fatality 


injury 


fatality 


injury 


intercept 




Q Q^j — 3.67 

-J.8(^_4 Q7 


-8.68_g:7g 


3.393_4 -[^4 


0^7 — 7.93 
-«.»'_9.99 


-3.94:3; J3 


-7 04^0.96 
'•y*_9.08 


-3.94:|;J^ 


wint 


_ ,09-129 


.iyz_ 256 






_ ,^„-.0641 
.loa_ 262 


_ 1.Q-.0641 
.j.o»_ 262 


_ -,^„-.0641 
.103_ 262 


-2i7i:°r 


jobend 


141-208 


141.208 
•-^*-^.0730 


142-209 


142-209 
■^^^^.0760 


- 


144.212 
•-^*^.076B 


- 


•i44:g^i5 


cons 


-.27o::°^32 


Qyri — .0532 
■^"J-.487 


- 97q--0B44 
■^'^-.496 


-■279::°^f 


cy QQ — ■ 393 

— ^-^^-5.07 


- 


Q — . 393 

— 07 


-.5981,202 


day 


-■7791^^^ 


.0654:i;«3 


•' '"-1.03 


.0784:131^ 


-.768:,^« 


- 


-768:5,« 


1 '?q.251 

-^■^^.0329 


snow 


-1-92I3'33 


-•3701:491 


10-. 861 


_ o74-.2,54 
496 


Q ~ ■ 265 
5J2 


-•388:;2f5 


-..388I:^«5 


Q Q — ■ 265 
-•3»»_.512 


dry 


.00 r 264 


•299:i« 


tr7a.887 
■3'».281 


-298:ig 


7151.02 
418 


.297:359 


71 5I-O2 
- ' l^°.418 


■297:359 


way4 


■308; 236 


.308:236 


ono. 376 
.ouo 231 


•303:ii? 


01 q. 433 
.J19 205 


01Q.433 

.c>iy 205 


- 


.308:111 


driver 




^■^^1.10 


io4.13 
■'•-■-''2.30 


l-18Hg 


104-14 
•^-l'^2.26 


i-27i;?i 


1 271-39 
l-^'l.l5 


1-0418^1 


MIO 


979-'>33 


■789:»« 


ICC. 433 
■-'-"'^-.0966 


.811:11;] 




-807;?6| 




■807:?tl 


moto 


2 -532:15 


2 -531: 35 


9 542.72 


9 542.72 

^•°^2.36 


2 552-73 
Z.OO2.37 


2.552;|3 


2-55|P 


2.55|:37 


vagc 


Qoi 9.0358 


QOI 9.0358 


noi 9.03,58 


n-^l 9.0358 


-0348:0411 


■0348;04ll 


-0348:0411 


.0249:0334 


-V27 










■(«io:il!;:ii„ 


(,., in.lK.ii 
-"■'-'-",UIJ2!)!) 


-2i3:i^H 


-2i-^;i25 


Ind 


•36i:i; 


•36i:i'i' 


.359:115 


-359:i[; 


.3C9.463 
.JDZ 263 


.ODZ 263 




.ouz 263 


X29 


6.08i:?f 
X 10~^ 


6.08i:?f 
X 10~^ 


6^3oi:i» 

X 10~^ 


6-3oi:io 

X 10~^ 




fi 949- 15 
''■^^3.30 

X 10~^ 




6 949 -15 
"•^*3.30 

X 10"^ 


priv 


— fi7q~'^'*2 

■0'y-.852 


-•679::|g 


.DMZ_ g4g 


-692I:g^ 


Q yc-- 1.73 
•J- '"'-6.55 


-3.6591 :|?4 


-3.7511;^ 


-3.659i:«;4 


X33 




•819i5|jl 


1.932;j^2 


-82515?^ 


2 4Q3.21 
^-^■^1.69 


•8O8I5OJ 




■80815"^ 


^34 


Q1 OA. 0202 


•00318:l»f™ 


ni on. 0200 


.00318:|»f™ 


(11 45-0215 
.Ui^O 00719 




0145-0215 
-'^1^''. 00719 


00692-00998 
.uuoyz 00396 


^35 




■286:233 




000.341 

.^»» 234 


^nK--211 
-.495_ 


.292:239 


-.495I:?iJ 


909.345 
.zaz 239 


driv 




007.440 
■'^°' .333 




.385.438.331 


■398:320 


QqQ.47B 
.CS90 320 




317-421 
' .209 








.00858 


.309 


.0695 


.293 


.0115 


.335 








909.428 

•^''^.140 








.436:652 


po and pi 






.607:^32 and .393:ii 


# free par. 


29 


29 


36 


averaged LL 




inncQ QQ-19046.91 
l»UOO.O»_jgQ6j 68 


-18952.63i;^835.03 


max(LL) 


-19023.62 (true) 


-19041.28 (observed) 


-18915.07 (observed) 


marginal LL 




-19065 07-19061.88 


-18994.00:JS:7i 


Good.-of-fit 




0.398 


0.601 


max(PSRF) 




1.00267 


1.00055 


MPSRF 




1.00310 


1.00073 


# observ. 


accidents = fatalities + injuries + PDOs: 32236 = 281 + 9947 + 22008 



Tabic 7.6 

Estimation results for multinomial logit models of severity outcomes of two- vehicle accidents on Indiana streets 



Variable 


ML-by-MLE 


ML-by-MCMC 


MSML 


state s = 


state s = 1 


fatality 


injury 


fatality 


injury 


fatality 


injury 


fatality 


injury 




in «~9-58 
-10.b_ii.6 


Qfi — 2.71 


-10.7_ii 7 


„ Qt:-2.79 
^•»O_3.10 


-13.1_ig 2 


nn — 2.87 
-3.UU_3 12 


1 1 -11-0 


Q nn — 2.87 
-3.00_3 12 






IOC-. 101 


10/ -.0999 


1O/I-.0999 




ion -.0939 
^•13" -.165 




1 on -.0939 
^■-■-■^"-.165 


'^^v 




ini 


^■•^''--1,24 


1(10 ■"''79 


^■'^•^•'-i,l« 




^■'^•^•'-i.l« 


-.U,)M)_ ^32 




.UOOU_ ggj^3 


-.U&5U_ (,<js3 


.nt--. 00559 
-.4»t>_ (,,jjg 


00559 
-.4»tl_ 0916 




— .ubt)y_ 121 








rvofll ~ 0188 

-.U8Ui_ j^42 


nsni ^ .0188 
-.USUi_ j^42 


ncQO-. 00109 


r,c:QO-. 00109 

-.u&y8_ 120 












14fi- 0465 
-.i4b_ 246 


-.i4b_ 246 


, . .-.0455 
-.i44_ 244 


1 ~-0455 
-•i44_ 244 




1OQ--0411 
-■13y_ 239 




ioq-.0411 
-.13y_ 239 




1 QQ-237 


1 QQ-237 

.lyy j^62 


1 QA-232 
■1^4.156 


1 Qzl-232 


1 n'^i-38 

l.uo 672 


1 OQ-226 
•1°°.150 


1 n'^1-38 

l.uo 672 


1 QQ .226 
•1°°.150 




-.282I:2f2 










rjQQ — .243 

324 


-.2721:188 


-.2721:188 


TioTi roaH 




-.6542,^9 


-697Iii«o 


-697Iii«o 


- 697-191 
■oy ' _1.26 


-.697Iii9i 


- 697- 191 
■oy '-1.26 


-.697Iii9i 


hllO 


.763-111 


.763:^5 


■802:?ii 


.802:863 


.801:?35 


-8oi:?i 


.801:835 


.801:835 


moto 




l-76l:i 


4.66t:i» 


1.751:11 


4-66t:lf 


1.751:11 


4-66t:l? 


i-75l:i 


voldg 






•428:o885 


(-,004.0639 

.UOZ'i 000866 




n49c;-0805 
.U4Z0 00511 






Ind 


.0769:1305 


.0769:1305 


.0778:0^53 


.0778:0^53 


.0803:0271 


.0803:0271 


.0803:0271 


.0803:0271 


X29 


.081lio*.0580 


.0284:0307 


■081:0676 


.0286:0309 


•0797:1033 


0290-0312 
.uzyu 0267 


•0797:i03g 


0290-0312 
.uzyu 0267 


priv 






_ 1:40-. 400 
■'^^'^-.689 


-543I:«o 


_ 500- 396 

•■^'^■-'-.685 


_ .00- 396 
•"J ''^-.685 


_ 5oq--396 


_ .oq-.396 
-"^■'^-.685 




3.i4i:i 


1 ^cil-73 
^•°>^1.37 


3 n7»-*^i 




1^°^1.30 


1 54I-8I 
1-^^1.30 


1 541-81 

1-°^1.30 


i.7o?:« 


X34 


.UIDZ QQ732 




n-if;n.0248 
■'^^O'J. 00714 




(1170-0268 
.Ul (y 00881 




r|i7q.0268 
.Ul (y 00881 




singTR 


7771.33 
•'".221 


.010_ 3gg 


-.758lj7^ 


— .310_ 3g2 


.you 300 


-.3061:235 




-.3061:235 


maxpass 


n'i9fi0615 
.UOZD 0437 


nc;9f:.0615 


.UO/iO 043g 


ncoaOeiS 
.uozo 0433 


.Ujyo 0292 


0*^08 -0501 

.UCjyO Q2g2 


.Ujyo 0292 


153-192 
.100 120 


mm 


cQi .926 
.001 236 


-23o::»? 


.582:237 


990-197 
260 


rroq.883 
•■^■^^.195 


-•26o::^j^ 


.539:883 


,oc;-.0500 
.±O.J_ 216 


slush 




.ZU41_ 300 




- 211-115 
•^^^-.307 




-•207I:iii 




-•207I:iii 



CO 
-J 



Table 7.6: (Continued) 



Variable 


ML-by-MLE 


ML-by-MCMC 


MSML 


state s = 


state s = 1 


fatality 


injury 


fatality 


injury 


fatality 


injury 


fatality 


injury 


driver 




■1 '''2; 0856 




•-^ ' ^.0859 


2.07^2i6 


164-237 
--'-""'-.0900 


2.07^^2?^ 




A27 




nice:-. 00346 
-.0165_ 0296 




.UiD.3_ Q2g3 








.UZUCi_ 0341 


nosig 




-.186_;223 




_ 104-158 
.iy4_ 




.ia4_ 230 




-•194I:18» 


smgo U V 




-.0860I:?ff 




-.08541:05™ 




-.0864I:°ff 




-.08641:0588 


oldvage 




0205-0236 
.UZUO oj^74 




020^ 0236 
.UZUO „j74 


.0205:0f35 


.0205:0235 


.0205:0235 


.UZUO 0175 


cigeOo 




- 521 - ^^-^S 




P.29-.349 
■"-'^^-.701 


— f^9fi — ■'^52 
•'-'^"-.706 


-.5261:382 


— f^9fi— -382 


— ^9fi— -352 
■^^'-'-.706 








.00107 


.221 


.00112 


.218 


.00091 


.232 








91 7.360 
•^-'-'.107 


i'i— 








po and pi 






.733:111 and .267:fii 


# free par. 


36 


36 


39 


averaged LL 




-64232.05:«St:i 


c/licro riT — f'4134.19 
-64152. 07_f;4i72. 22 


max(LL) 


-64226.29 (true) 


-64217.50 (observed) 


-64113.04 (observed) 


marginal LL 




-64245 77-64241.79 
DIZIO./ '_fi4247.82 


-64191 23-64180.82 

"^-^='-^■^'^-64193, 80 


(iood.-of-Iil 




0.77:-! 


0.781 


max(PSRF) 




1.00092 


1.00569 


MPSRF 




1.00152 


1.00658 


# observ. 


accidents = fatalities + injuries + PDOs; 125336 = 138 + 27727 + 97471 



Table 7.7 

Explanations and summary statistics for variables and parameters listed in Tables [7111 - 17.61 and in Tables IXn - |A.4l 



Variable 


Description 


Mean 


Std'' 


Min'' 


Median 


Max'' 


agcO 


Age of the driver at fault is less than 18 years old (dummy) 


.0846 


.278 








1.00 


agoOo 


Ago of the oldest driver involved into the accident is less than 18 years old (dummy) 


.0103 


.101 








1.00 


cons 


Construction at the accident location (dummy) 


.0272 


.163 








1.00 


curve 


Roadway is at curve (dummy) 


.0459 


.209 








1.00 


dark 


Dark time with no street lights (dummy) 


.0439 


.205 








1.00 


darklamp 


Dark and street lights on (dummy) 


.130 


.337 








1.00 


day 


Daylight (dummy) 


.784 


.412 





1.00 


1.00 


dayt 


Day hours: 9:00 to 17:00 (dummy) 


.577 


.495 





1.00 


1.00 


driv 


Roadway median is drivable (dummy) 


.415 


.493 








1.00 


driver 


Primary cause of the accident is driver-related (dummy) 


.964 


.185 





1.00 


1.00 


dry 


Roadway surface is dry (dummy) 


.739 


.439 





1.00 


1.00 


env 


Primary cause of the accident is environment-related (dummy) 


.0255 


.158 








1.00 


lillO 


Help arrived in 10 minutes or less after the crash (dummy) 


.637 


.481 





1.00 


1.00 


hl20 


Help arrived in 20 minutes or less after the crash (dummy) 


.834 


.372 





1.00 


1.00 


Ind 


License state of the vehicle at fault is Indiana (dummy) 


.907 


.290 





1.00 


1.00 


liglit 


Daylight or street lights are lit up if dark (dummy) 


.914 


.281 





1.00 


1.00 


maxpass 


The largest number of occupants in all vehicles involved 


1.88 


1.77 







70.0 


mm 


Two male drivers are involved, if a 2-vehicle accident (dummy) 


.308 


.461 








1.00 


morn 


Morning hours: 5:00 to 9:00 (dummy) 


.131 


.337 








1.00 


moto 


The vehicle at fault is a motorcycle (dummy) 


.00348 


.0589 








1.00 


nigh 


Late night hours: 1:00 to 5:00 (dummy) 


.0148 


.121 








1.00 


nocons 


No construction at the accident location (dummy) 


.973 


.163 





1.00 


1.00 


nojun 


No roadway junction at the accident location (dummy) 


.448 


.497 





1.00 


1.00 


nonroad 


Non-roadway crash (parking lot, etc.) (dummy) 


.00518 


.0718 








1.00 



Table 7.7: (Continued) 



Variable 


Description 


Mean 


Std'' 


Min'' 


Median 


Max'' 


nosig 


No any traffic control device for the vehicle at fault (dummy) 


.233 


.423 








1.00 


olddrv 


The driver at fault is older than the other driver, if a 2-vehicle accident (dummy) 


47.3 


16.5 


15.0 




99.0 


oldvage 


Age of the oldest vehicle involved (in years) 


10.2 


5.07 


-1.00 




41.0 


othUS 


License state of the vehicle at fault is a U.S. state except Indiana and its neighboring 
states (IL, KY, OH, MI) (dummy) 


.0272 


.148 








1.00 


precip 


Precipitation: rain/freezing rain/snow/sleet/hail (dummy) 


.172 


.377 








1.00 


priv 


Road traveled by the vehicle at fault is a private drive (dummy) 


.0289 


.168 








1.00 


r21 


Roadway traveled by the vehicle at fault is two-lane and one-way (dummy) 


.0347 


.183 








1.00 


rmd2 


Roadway traveled by the vehicle at fault is multi-lane and divided two-way (dummy) 


.230 


.421 








1.00 


SingSUV 


One of the two vehicles involved is a pickup OR a van OR a sport utility vehicle, 

if a. 2-\'C'liit*ie a.tx'iciciil (ciuiiiiiu'j 


.111) 


.197 








1.00 


singTR 


One of the two vehicles is a truck OR a tractor, if a 2-vehicle accident (dummy) 


.0688 


.253 








1.00 


slush 


Roadway surface is covered by snow/slush (dummy) 


.0400 


.196 








1.00 


snow 


Snowing weather (dummy) 


.0414 


.199 








1.00 


str 


Roadway is straight (dummy) 


.949 


.220 





1.00 


1.00 


sum 


Summer season (dummy) 


.243 


.429 








1.00 


sund 


Sunday (dummy) 


.0784 


.269 








1.00 


thday 


Thursday (dummy) 


.157 


.364 








1.00 


vage 


Age of the vehicle at fault (in years) 


7.91 


5.31 


-1.00 




41.0 


voldg 


The vehicle at fault is more than 7 years old (dummy) 


.489 


.500 








1.00 


voldo 


Age of the oldest vehicle involved is more than 7 years (dummy) 


.688 


.463 





1.00 


1.00 


wall 


Road median is a wall (dummy) 


.0528 


.224 








1.00 


way4 


Accident location is at a 4-way intersection (dummy) 


.371 


.483 








1.00 


wday 


Weekday (Monday through Friday) (dummy) 


.800 


.400 





1.00 


1.00 


wint 


Winter season (dummy) 


.250 


.433 








1.00 



o 
o 



Table 7.7: (Continued) 



Variable 


Description 


Mean 


Std'' 


Min'' 


Median 


Max"" 


Xi2 


Roadway type (dummy: 1 if urban, if rural) 


.829 


.377 





1.00 


1.00 


X27 


Number of occupants in the vehicle at fault 


1.45 


1.18 







70.0 


X2g 


Speed limit (used if known and the same for all vehicles involved) 


36.7 


9.86 


5.00 




75.0 


^33 


At least one of the vehicles involved was on fire (dummy) 


.00505 


.0709 








1.00 




Age of the driver at fault (in years) 


37.0 


9.86 


3.00 




99.0 


X35 


Gender of the driver at fault (dummy: 1 if female, if male) 


.449 


.497 








1.00 




Probability of i^^ severity outcome averaged over all values of explanatory variables Xt_„ 


_ 


_ 


_ 


_ 


_ 




Markov transition probability of jump — > 1, as time t increases to t + 1 














Markov transition probability of jump 1 — ► 0, as time t increases to t + 1 












po and pi 


Unconditional probabilities of states and 1 












# free par. 


Total number of free model parameters ()3-s) 












averaged LL 


Posterior average of the log-likelihood (LL) 












max(LL) 


True maximum value of log-likelihood (LL) for MLE; maximum observed value of LL 
for Bayesian-MCMC 












marginal LL 


Logarithm of marginal likelihood of data (ln[/(Y A^)]) 












Good.-of-fit 


Goodness-of-fit p-value, refer to equation l|4.5|l 












max(PSRF) 


Maximum of the potential scale reduction factors ^ 












MPSRF 


Multivariate potential scale reduction factor (MPSRF) ^ 












# observ. 


number of observations of accident severity outcomes available in the data sample 













^ Standard deviation, minimum and maximum of a variable. 

^ PSRF/MPSRF are calculated separately/jointly for all continuous model parameters. PSRF and MPSRF are close to 1 for converged MCMC chains. 
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From Tables I7.1H7.6I we find that in all cases when the two states and Markov 
switching multinomial logit (MSML) models exist, these models are strongly favored 
by the empirical data over the corresponding standard multinomial logit (ML) models. 
Indeed, for example, from lines "marginal LL" in Tables I7.1H7.6I we see that the 
MSML models provide considerable, ranging from 40.49 to 206.41, improvements of 
the logarithm of the marginal likelihood of the data as compared to the corresponding 
ML models^ Thus, from equation (14. 3 p we find that, given the accident severity data, 
the posterior probabilities of the MSML models are larger than the probabilities of 
the corresponding ML models by factors ranging from e^°'^^ to e^''^'^^. Note that we 
use equation (14. 2 p for calculation of the values and the 95% confidence intervals of 
the logarithms of the marginal likelihoods. The confidence intervals are found by 
bootstrap simulations (see footnote [7] on page [621) . 

Note that a classical statistics approach for model comparison, based on the max- 
imum likelihood estimation (MLE), also favors the MSML models over the standard 
ML models. For example, refer to line "max(LL)" in Table mi given for the case of 1- 
vehicle accidents on interstate highways. The MLE gave the maximum log-likelihood 
value —8465.79 for the standard ML model. The maximum log-likelihood value ob- 
served during our MCMC simulations for the MSML model is equal to —8358.97. 
An imaginary MLE, at its convergence, would give a MSML log-likelihood value that 
would be even larger than this observed value. Therefore, if estimated by the MLE, 
the MSML model would provide large, at least 106.82 improvement in the maximum 
log-likelihood value over the corresponding ML model. This improvement would come 
with only modest increase in the number of free continuous model parameters (/?-s) 
that enter the likelihood function (refer to Table [7?T] under free par."). Simi- 
lar arguments hold for comparison of MSML and ML models estimated for other 
roadway-class-accident-type combinations where two states of roadway safety exist 
(see Tables OHESD- 

"^In addition, we find that DIG (deviance information criterion) favors the MSML models over the 
corresponding ML models by DIG value improvement ranging from 168.33 to 450.52. However, we 
prefer to rely on the Bayes factor approach instead of the DIG (see footnote [2] on page [3T|) . 
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Figure 7.1. Weekly posterior probabilities P{st = 1|Y) for the MSML 
models estimated for severity of 1- vehicle accidents on interstate highways 
(top plot), US routes (middle plot) and state routes (bottom plot). 



To evaluate the goodness-of-fit for a model, we use the posterior (or MLE) es- 
timates of all continuous model parameters (/3-s, a, Po^i, Pi^o) and generate lO'' 
artificial data sets under the hypothesis that the model is true (see footnote [T7] on 
page[H3]). We find the distribution of x^, given by equation 04.51) . and calculate the 
goodness-of-fit p- value for the observed value of x^- The resulting p- values for our 
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Jan-03 Jul-03 Jan-04 Jul-04 Jan-05 Jul-05 Jan-06 Jul-06 

Date 



Figure 7.2. Weekly posterior probabilities P{st = 1|Y) for the MSML 
models estimated for severity of l-vehicle accidents occurring on county 
roads (top plot), streets (middle plot) and for 2- vehicle accidents occurring 
on streets (bottom plot). 

multinomial logit models are given in Tables I7.1H7.6I These p-values are around 
20-80%. Therefore, all models fit the data well. 

Now, refer to Table 17.81 The first six rows of this table list time-correlation 
coefficients between posterior probabilities P{st = 1|Y) for the six MSML models that 
exist and are estimated for six roadway-class- accident-type combinations (l-vehicle 
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accidents on interstate highways, US routes, state routes, county roads, streets, and 
2- vehicle accidents on streets) We see that the states for 1- vehicle accidents on all 
high-speed roads (interstate highways, US routes, state routes and county roads) are 
correlated with each other. The values of the corresponding correlation coefficients 
are positive and range from 0.263 to 0.688 (see Table 17781) . This result suggests an 
existence of common (unobservable) factors that can cause switching between states 
of roadway safety for 1-vehicle accidents on all high-speed roads. 

The remaining rows of Table 17.81 show correlation coefficients between posterior 
probabilities P{st = 1|Y) and weather-condition variables. These correlations were 
found by using daily and hourly historical weather data in Indiana, available at the 
Indiana State Climate Office at Purdue University (www.agry.purdue.edu/climate). 
For these correlations, the precipitation and snowfall amounts are daily amounts in 
inches averaged over the week and across Indiana weather observation stations (see 
footnote [19] on page [85]) . The temperature variable is the mean daily air temperature 
{°F) averaged over the week and across the weather stations. The wind gust vari- 
able is the maximal instantaneous wind speed (mph) measured during the 10-minute 
period just prior to the observational time. Wind gusts are measured every hour 
and averaged over the week and across the weather stations. The effect of fog/frost 
is captured by a dummy variable that is equal to one if and only if the difference 
between air and dewpoint temperatures does not exceed 5°F (in this case frost can 
form if the dewpoint is below the freezing point 32° F, and fog can form otherwise). 
The fog/frost dummies are calculated for every hour and are averaged over the week 
and across the weather stations. Finally, visibility distance variable is the harmonic 
mean of hourly visibility distances, which are measured in miles every hour and are 
averaged over the week and across the weather stations (see footnote [20] on page [861) . 

From the results given in Table [778] we find that for 1-vehicle accidents on all high- 
speed roads (interstate highways, US routes, state routes and county roads), the less 
frequent state = 1 is positively correlated with extreme temperatures (low during 



'See footnote [13] on page [77] for details on computation of correlation coefficients. 
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Table 7.8 

Correlations of the posterior probabilities P{st = 1|Y) with each other 
and with weather-condition variables (for the MSML models of accident 
severities) 





1-vehicle, 


1-vehicle, 


1-vehicle, 


1-vehicle, 


1-vehicle, 


2-vehicle, 




interstates 


US routes 


state routes 


county roads 


streets 


streets 


l-\'C'iiit*ic. iiilc'isla.lc's 


1 


0. 118 


0.29:-! 


().()()() 


-o.oi:i 


-0.173 


1-vehicle, US routes 


0.418 


1 
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-0.070 


-0.155 
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0.293 
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1 


0.409 


-0.047 


-0.035 


1-vehicle, county roads 
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1 


-0.022 


-0.051 


1-vehicle, streets 


-0.013 
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-0.022 


1 


0.115 
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-0.173 


-0.155 


-0.035 


-0.051 


0.115 


1 
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Precipitation (inch) 
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Temperature {°F) 
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-0.439 
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> 0.0 (dummy) 
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0.164 
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0.167 
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Precipitation (inch) 
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Temperature {°F) 
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-0.124 
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0.077 
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winter and high during summer), rain precipitations and snowfalls, strong wind gusts, 
fogs and frosts, low visibility distances. It is reasonable to expect that roadway safety 
is different during bad weather as compared to better weather, resulting in the two- 
state nature of roadway safety. 

The results of Table 17^ suggest that Markov switching for road safety on streets is 
very different from switching on all other roadway classes. In particular, the states of 
roadway safety on streets exhibit low correlation with states on other roads. In addi- 
tion, only streets exhibit Markov switching in the case of 2-vehicle accidents. Finally, 
states of roadway safety on streets show little correlation with weather conditions. A 
possible explanation of these differences is that streets are mostly located in urban 
areas and they have traffic moving at speeds lower that those on other roads. 

Next, we consider the estimation results for the stationary unconditional proba- 
bilities pQ and pi of states St = and St = I for MSML models [see equations fl3.16p ]. 
These transition probabilities are listed in lines "po and pi" of Tables I7.1H7.6I We 
find that the ratio pi/po is approximately equal to 0.46, 0.13, 0.74, 0.25, 0.65 and 
0.36 in the cases of 1- vehicle accidents on interstate highways, US routes, state routes, 
county roads, streets, and 2-vehicle accidents on streets respectively. Thus, for some 
roadway-class-accident-type combinations (for example, 1-vehicle accidents on US 
routes) the less frequent state St = 1 is quite rare, while for other combinations (for 
example, 1-vehicle accidents on state routes) state = 1 is only slightly less frequent 
than state St = 0. 

Finally, we set model coefficients /3(o) and /3(-^-) to their posterior means, calcu- 
late the probabilities of fatality and injury outcomes in states and 1 by using 
equation (13.141) . and average these probabilities over all values of the explanatory 
variables Xj „ observed in the data sample. We compare these probabilities across 
the two states of roadway safety, St = and St = 1, for MSML models [refer to lines 
"(P/*2)x" Tables [7n - l7. 6] . We find that in many cases these averaged probabilities 
of fatality and injury outcomes do not differ very significantly across the two states 
of roadway safety (the only significant differences are for fatality probabilities in the 
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cases of 1-vehicle accidents on US routes, county roads and streets). This means that 
in many cases states St — and St — 1 are approximately equally dangerous as far 

as accident severity is concerned. We discuss this result in the next chapter (which 
includes a discussion of all our results) . 
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CHAPTER 8. SUMMARY AND CONCLUSIONS 

In this final chapter we give our major conclusions for the two-state Markov switching 
models estimated for annual accident frequencies, weekly accident frequencies, and 
for accident severities. 

• Our conclusions for the Markov switching models of annual accident frequencies, 
specified in Section 13.41 and estimated in Section 16.11 are as follows. First, 
these models provide a far superior statistical fit for accident frequencies as 
compared to the standard zero-inflated models. Second, the Markov switching 
models explicitly consider transitions between the zero-accident state and the 
unsafe state over time, and permit a direct empirical estimation of what states 
roadway segments are in at different time periods. In particular, we found 
evidence that some roadway segments changed their states over time (see the 
bottom-right plot in Figure [6^ . Third, note that the Markov switching models 
avoid a theoretically implausible assumption that some roadway segments are 
always safe because, in these models, any segment has a non-zero probability 
of being in the unsafe state. Indeed, the long-term unconditional mean of 
the accident rate for the n*^ roadway segment is equal to Pi'\\t,n)t-, where 
Pi"'' = Po"^i/(Po'^i + Pi"!o) is stationary probability of being in the unsafe 
state St^n = 1 and {Xt,n)t is the time average of the accident rate in the unsafe 
state [refer to equations (13.71) and (13.161) ]. This long-term mean is always above 
zero (see the bottom plot in Figure . even for segments that seem to be in 
the zero-accident state over the whole observed five-year time interval of our 
empirical data. Finally, we conclude that two-state Markov switching count 
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data models are likely to be a better alternative to zero-inflated models, in 
order to account for excess of zeros observed in accident frequency data. 

Our conclusions for the Markov switching models of weekly accident frequen- 
cies, specifled in Section 13.51 and estimated in Section 16. 2[ follows. Our 
empirical flnding that two states exist and that these states are correlated with 
weather conditions has important implications. For example, multiple states of 
roadway safety can potentially exist due to slow and/or inadequate adjustment 
by drivers (and possibly by roadway maintenance services) to adverse conditions 
and other unpredictable, unidentified, and/or unobservable variables that influ- 
ence roadway safety. All these variables are likely to interact and change over 
time, resulting in transitions from one state to another. As discussed earlier, the 
empirical findings show that the less frequent state is significantly less safe than 
the other, more frequent state. The estimation results of the full MSNB/MSP 
models show that explanatory variables ^t,n exert different influences on road- 
way safety in different states as indicated by the fact that some of the parameter 
estimates for the two states of the full MSNB/MSP models are signiflcantly dif- 
ferent. Thus, the states not only differ by average accident frequencies, but also 
differ in the magnitude and/or direction of the effects that various variables 
exert on accident frequencies. This again underscores the importance of the 
two-state approach^ 

Our conclusions for the Markov switching models of accident severities, specified 
in Section 13.61 and estimated in Chapter [71 dbTQ clS follows. We found that 
two states of roadway safety and Markov switching multinomial logit (MSML) 
models exist for severity of 1-vehicle accidents occurring on high-speed roads 
(interstate highways, US routes, state routes, county roads), but not for 2- 
vehicle accidents on these roads. One of possible explanations of this result 



^ One might also consider a threshold model in which the state value is a function o f explanatory vari- 



ables [similar to threshold autoregressive models used in econometrics [TsajJ,[2002|]. This interesting 
possibility is beyond the scope of this study. 



Ill 



is that 1- and 2-vehicle accidents may differ in their nature. For example, 
on one hand, severity of 1-vehicle accidents may frequently be determined by 
driver-related factors (speeding, falling a sleep, driving under the influence, 
etc). Drivers' behavior might exhibit a two-state pattern. In particular, drivers 
might be overconfident and/or have difficulties in adjustments to bad weather 
conditions. On the other hand, severity of a 2-vehicle accident might crucially 
depend on the actual physics involved in the collision between the two cars (for 
example, head-on and side impacts are more dangerous than rear-end collisions). 
As far as slow-speed streets are concerned, in this case both 1- and 2-vehicle 
accidents exhibit two-state nature for their severity. Further studies are needed 
to understand these results. In this study, the important result is that in all 
cases when two states of roadway safety exist, the two-state MSML models 
provide a superior statistical fit for accident severity outcomes as compared to 
the standard ML models. 

We found that in many cases states St — and St — 1 are approximately equally 
dangerous as far as accident severity is concerned. This result holds despite the 
fact that state St = 1 is correlated with adverse weather conditions. A likely 
and simple explanation of this finding is that during bad weather both num- 
ber of serious accidents (fatalities and injuries) and number of minor accidents 
(PDOs) increase, so that their relative fraction stays approximately constant. In 
addition, most drivers are rational and they are likely to take some precautions 
while driving during bad weather. From the results of modeling annual accident 
frequencies, we know that the total number of accidents significantly increases 
during adverse weather conditions. Thus, driver's precautions are probably not 
sufficient to avoid increases in accident rates during bad weather. 

We can speculate that one of the major causes of the existence of different states 
of roadway safety can be slow and inadequate adjustment by some drivers to sudden 
worsening of weather and roadway conditions (such as snow or ice on a roadway). 
Of course, apart from weather conditions, there can be additional unpredictable and 
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unidentified factors tliat influence road safety. All these factors are likely to interact 
and change in time, resulting in unobserved heterogeneity in accident data. Markov 
switching between states of roadway safety intends to account for these factors and 
for the resulting unobserved heterogeneityo Examples of other statistical models that 



intend to account for unobserved heterogeneity, include finite mixture mode- 
mixed ) models, and randon i effects models 



dom parameters 



Washington et al. 



Shankar et al 



s , ran- 



2003, 



Park and Lord 



1998 . 



20081 ]. A theoretical advantage of Markov 



switching models over other models is that the former allows for an explicit identifi- 
cation of the states of roadway safety at different time periods. Another advantage 
of Markov switching models is that they explicitly consider how various explanatory 
variables exert different influences on road safety in different states. For example, 
in the case of the MSNB and MSP models of accident frequencies estimated in this 
study, the states differ not only by the values of the average accident frequency (A), 
but also by the values of the model coefficients (/?-s) in the two states. 

As far as practical application of Markov switching models for prediction of aver- 
aged accident rates is concerned, this prediction depends on whether it is conditional 
or unconditional. For probabilities conditioned on the previous state, one uses the 
transition probabilitieso For all unconditional expectations and long-term predic- 
tions, one uses unconditional probabilities (pg"''' and Pi"^), given by equation (13.161) . 
In particular, the long-term probability of being in a state is equal to the unconditional 
probability of this state. Please note that, even if the current state is known (zero or 
one), then in a long run, all expectations converge to the unconditional expectations 
exponentially fast (this is a property of Markov processes). Because researchers and 
practitioners are usually interested in a long-term improvement of safety, using the 
unconditional probabilities is more appropriate for predictions and decision making. 



^The Markov property of the switching serves as a reasonable approximation, which helps to sim- 
plify our analysis. For example, the Markov property holds reasonably well for changes of weather 
conditions in time. 

■^For example, if the previous state was zero, Sf-i.„ — 0, then the probabilities of the current state 
St^n being zero and one are equal to the transition probabilities Pq!!^q and p'^i respectively, refer to 
equation p.lSp . 
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A determination of the roadway safety state value (zero or one) during a specific 
time period t is complicated by the unobservability of the state variable. As a result, 
we rely on Bayesian inference in this case - we use an accident data, estimate a Markov 
switching model for this data, and find the posterior probabilities for the state values 
at time t. These posterior probabilities should be used for inference about the state 
values. 

In terms of future work on Markov switching models for accident frequencies and 

severities, additional empirical studies (for other accident data samples) and multi- 
state models (with more than two states of roadway safety) are two areas that would 
further demonstrate the potential of the approach. 
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Table A.l 

Estimation results for multinomial logit models of severity outcomes of 
two-vehicle accidents on Indiana interstate highways 
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Table A. 2 

Estimation results for multinomial logit models of severity outcomes of 
two- vehicle accidents on Indiana US routes 
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-■4221:^!^? 


1 91 - .0619 
~-l^l-.179 


Xl2 


1 QC— .955 
l.oo_j 75 


— .313_ 3gg 




— 'Hon— -2-1* 

.ozu_ 392 


dark 


.546:?ii 


115-229 
■^-'■'^-.00220 


^40.926 


00229 


snow 




- 259- "''0''* 

.ZO.J_ 


r,f.r,-.0952 

.ZDZ_ 43J 


9f.9~.O952 
.ZDZ_ 43J 


driv 


.UDUU_ QQ240 


.UDUU_ QQ24g 


■0556- j;?0157 


.0556i11qqj^57 


nojun 


or>9.582 


•^-'-^-.269 


nono.583 




driver 


•426:iS 


•^^",280 


•428::l:i 


■428:ii 


hllO 




)8() 


.r)(ii:gg 




moto 


3 -983: Si 


1 aa2.24 
1-°°1.61 


3-97|:i? 


1.882:28 


vage 






•0482:«™t 






0740.0999 


noqi .0268 


•0757:ioii 


n9Qo.0270 

.UZOO gjgg 


priv 


_1 1Q--540 


_1 1Q--540 
73 


1 1O-.607 
-1-10_1.81 


1 10-. 607 
-1-18_1.81 


-'^33 


2-981:32 




9 Q73.62 
^■y'2.28 


1-^^1.03 


singTR 






1 151-44 
-■-■^^.843 




maxpass 


-0776:0|?1 


.0776:{5i?l 


■0784:0111 


-0784:0i|i 


olddrv 


-0198:S 


■0230:;|?«f 


•0199:0??fi 


•00481:«jitt 


mm 


01 p. 598 
••^1°.0343 


•00480:00|i0 


•391 .602 
••^^1.0417 


— .230_ 289 


oldvage 




noi— .175 

.zo'i_ 292 




.0230:gfff 


{pI%)x 






-00759 


.255 














Po and pi 






# free par. 


32 


32 


averaged LL 




— Ifi5'^5 4c;-16628.62 


max(LL) 


-16527-94 (true) 


16522.89 (observed) 


marginal LL 




-16549.59i|i|t:iO 


Good.-of-fit 




0.372 


max(PSRF) 




1.00275 


MPSRF 




1.00358 


# observ. 


accid.=fatal-+inj.+PDO: 28259 = 222 + 7285 + 21022 
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Table A. 3 

Estimation results for multinomial logit models of severity outcomes of 
two- vehicle accidents on Indiana state routes 



Variable 


ML-by-MLE 


ML-by-MCMC 


fatality 


injury 


fatality 


injury 


intercept 


-IS 


Q ^[-—3.37 


-i3.2:;ij 


'J-'O_4.03 


wint 


-.0668I:?»^9» 


-. 06681 


-. 06691 


-.06691;°^^**^ 


wday 


-,00-. 0737 
■ ^^■^■^-.192 


-i33::Sf 


1 qo-.0727 
•J^'^^-.191 


1 00 -.0727 
■-'-•-'^-.191 


Xl2 




-25ii:^S 


-796Ii*?I 


-•2621 


dark 




.248:fii 


i.ori^il 


•248;?ii 


wall 






-2.56l5^°| 




nojun 


.385:^12 


— 1 7n~l'^l 


j^42 


-.172::12-J 


curve 


^•^^.715 


004.323 


i-ooi„^l 


•239:?i,^ 


driver 


1 nyl-Bs 

^^•^'.450 


400.542 
■^^^.301 


-'-■-'--'-.,521 


.4i8:i« 


hI20 


1 21 1"-* 


70c:. 810 


1 22l-'^" 

-'-•^^,780 


00c:. 981 
.««b 7S9 








— '-2.:u 




X29 


.0942:ii|4 


r,24fi.0277 


.0950;Jy|9 


•0249:°i«S^ 


priv 


-.856Ii^JS^ 


-.856 1 1^^^ 


-88IIS 


-ssilil,^ 


Xsi 


3.ioi:i 




3.io^:«t 


^•^'.950 


^35 


QQn.739 
.oou Q206 




•384;o324 




singTR 


1 

I.UUq 


- 114-°215 


i-ooifil 


T TO-. 0224 
206 


voldo 


.zoa 201 


ore. 309 
.zoa 201 


.ij* 200 


254-30* 
.za* 200 


maxpass 




n.;of5.0683 


Qt;44.0693 
.UO'1^1 0398 


-'J ''^^.0398 


olddrv 


noio.0284 


.UU40U 003m 


no-i 0.0284 
-'J^-'-^.014() 


n 4t;n-oo595 

.U.40U 00306 


mm 


.625;288 


'-.229 


coo. 975 
.uoo 305 


- 177--i2'i 

•^' '-.230 


nocons 




•280;f33 




■280:fi 


driver 




•^'^^.lee 




.460:^5 








.00843 


.257 














po and pi 






# free par. 


35 


35 


averaged LL 




oinaa -ji -21081.09 
-21U»».3i_2io97.38 


max(LL) 


-21096.20 (true) 


-21074.01 (observed) 


marginal LL 




-21103.71l2J09^;88 


Good.-of-fit 




0.635 


max(PSRF) 




1.00141 


MPSRF 




1.00176 


# observ. 


accid.=fatal.+inj.+PDO: 36136 = 311 + 9276 + 26549 
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Table A.4 

Estimation results for multinomial logit models of severity outcomes of 
two- vehicle accidents on Indiana county roads 



Variable 


ML-by-MLE 


ML-by-MCMC 


fatality 


injury 


fatality 


injury 


intercept 


-10.6I?,« 


— 3.50_3 72 


-in 7"'-'-'*i 
J^'J-'-ii.g 


c,o-3.37 

0.00_3 gg 


wint 


-.145I:°If 


.1410_ 214 


_ 146--°'^'^4 

.14:0_ 216 


.141D_ 216 


sund 


100.290 
■^^''.0945 


109-290 
•-'■='^.0945 


ion. 287 
■^^".0927 


ion. 287 
■^^".0927 


morn 


-.1081:?^^ 




■^"^-.181 


•^"^-.181 


X\2 


-1.48l2«I 


-i6o::°ir 




-1 65"-"*41 


darklamp 


_ -,07- 0239 

•^^'-.371 


•-'-^' -.371 


-.204I:!«42 


-.204::!«42 


way4 


.249:342 


.249:342 


.279:342 


.279:342 


driver 


.247:f^o 


.247:f^g 


.258:fi? 


.258:fi| 


hl20 


1 ';82ii 


Q14.993 




•9571304 


moto 


4 QA:'^-^'^ 




4.044:«I 


9 91 2.61 
Z.Zlj^ g2 


^29 


•0813:J«i5 


■0287:||ii!> 


.0820:i«i7 


n9Qn^0324 
■'J^yu.0257 




^•°^2,()(i 


l-18^7il! 


9 773.51 


-1 17I.56 

1.1 ( ygy 


slugs r\' 






■i-i:ES! 




oldvage 


■0390:8fi? 


noi c.()2fi9 
.UZIO „jg2 


.0387:8fli 


noi 7.0270 
•"^■^'.0163 


ageO 




1 49.230 
■-'-^^.0534 




•-'-^■-'.0552 


singTR 




_ 174- ■1)454 




-, VO-.0461 
— •-'^''^-.302 


maxpass 




ni 7c. 0286 
■"■"■'".00670 




017Q.0288 
■"■'■"'.00685 


ageOo 




C7C: — -335 
.0(0_ 




.OOd_ g2g 


mm 




— .258_ 322 




— .258_ 322 








.00662 


.247 


P0-.1 












po and Pi 






# free par. 


26 


26 


averaged LL 




-14423.80i;i«I:?^ 


max(LL) 


-14411.12 (true) 


-14412.78 (observed) 


marginal LL 




-144S4 7Q-14431.73 


Good.-of-fit 




0.370 


max(PSRF) 




1.00141 


MPSRF 




1.00225 


# observ. 


accid.=fatal.+inj.+PDO: 25597 = 173 + 6315 + 19109 



VITA 



122 



VITA 
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